Marine  Physical  Laboratory 


VLF  Source  Localization  with  a  Freely 
Drifting  Sensor  Array 

George  Chia-Jen  Chen 


SIO  Reference  92-18 


-OTIC 

tflELECTEf% 

OT  D 


MPL-U-45/92 
September  1992 


Approved  for  public  release;  Distribution  unlimited. 


University  of  California,  San  Diego 
Scripps  Institution  of  Oceanography 


DEFENSE  TECHNICAL  INFORMATION  CENTER 


6 


if  & 


nil  >w*ww 

9313593 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

PuMc  raportng  txsden  tor  this  eoHeetlon  of  Mormetion  I*  estlmeted  to  eversge  1  hour  par  rmtxnu,  Including  the  time  tor  reviewing  instructions,  sesrchmg  existing  deU  toixces. 
gathering  and  maintaining  the  data  need  ed,  and  completing  and  reviewing  the  cottectton  ot  information.  Send  comment*  regarding  thi*  burden  eattmate  or  any  other  aspect  ot 
this  eotiectlon  ot  mtormation,  including  suggestion*  lor  reducing  this  burden,  to  Washington  Headquarters  Services.  Directors!#  lor  intormstion  Operation*  and  Reports .  1S15  dederson 
□avia  Hkihwav.  Suite  1204.  Ahinqton.  VA  22202-4302  and  lo  the  OUce  ol  Manaoement  and  Sudoet,  Paperwork  Reduction  Protect  10704-0188).  Washirtqtor,  DC  20503. 

1.  Agency  Use  Only  (Leave  Blank).  2.  Report  Date.  3.  Report  Type  and  Dates  Covered. 

September  1992  Thesis 

4.  Title  and  Subtitle. 

VLF  Source  Localization  with  a  Freely  Drifting  Acoustic  Sensor 

Array 

5.  Funding  Numbers. 

N00014-88-K-2040 

Project  No. 

Task  No. 

6.  Author(s). 

George  Chla-Jen  Chen 

7.  Performing  Monitoring  Agency  Names(s)  and  Address(es). 

University  of  California,  San  Diego 

Marine  Physical  Laboratory 

Scripps  Institution  of  Oceanography 

San  Diego,  California  92152 

8.  Performing  organization 

Report  Number. 

MPL-U-45/92 

SIO  Reference  92-18 

9.  Sponsoring/Monitoring  Agency  Name(s)  and  Addresses). 

Naval  Research  Laboratory 

Atten:  Code  5160 

4555  Overlook  Avenue,  S.W. 

Washington,  D.C.  20375-5000 

10.  ^ons^ri|ijTgonitoring  Agency 

11.  Supplementary  Notes. 

12a.  Distribution/Availability  Statement. 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  Distribution  Code. 

13.  Abstract  (Maximum  200  words). 

Source  localization  and  tracking  capability  of  the  freely  drifting  Swallow  float  volumetric  array  is  dem¬ 
onstrated  with  matched-field  processing  (MFP)  technique  using  data  collected  during  the  1989  Swal¬ 
low  float  experiment  conducted  in  the  Northwest  Pacific.  Marine  Physical  Laboratory’s  set  of  nine 
freely  drifting,  infrasonic  sensors,  capable  of  recording  ocean  ambient  noise  in  the  1-  to  25-Hz  range, 
was  deployed  to  span  the  water  column  of  4100-m  depth,  with  horizontal  aperture  on  the  order  of  5  km. 
Even  though  the  floats  were  freely  drifting,  their  positions  were  determined  with  the  8-kHz  internal 
navigation  system  and  a  postprocessing  least-squares-based  localization  procedure. 

14.  Subject  Terms. 

Swallow  float,  matched-field  processing,  infrasonic  sensor,  vlf  source  localization 

16.  Price  Code. 

17.  Security  Classification  18.  Security  Classification  19.  Security  Classification 
or  Report  of  This  Page.  of  Abstract. 

Unclassified  Unclassified  Unclacsi.'ied 

20.  Limitation  of  Abstract. 

None 

NSN  754001-280-5500 


Standard  Form  298  (Rev.  2-89) 
Prescribed  by  ANSI  Sid.  239-18  298-102 


UNIVERSITY  OF  CALIFORNIA,  SAN  DIEGO 


VLF  Source  Localization  with  a  Freely  Drifting  Acoustic  Sensor  Array 


A  dissertation  submitted  in  partial  satisfaction  of  the 


requirements  of  the  degree  of  Doctor  of  Philosophy 


in  Electrical  Engineering  (Applied  Ocean  Science) 


by 


George  Chia-Jen  Chen 


Committee  in  Charge: 

Professor  William  S.  Hodgkiss,  Jr.,  Chair 
Professor  LeRoy  M.  Dorman 
Professor  John  A.  Hildebrand 
Professor  Laurence  B.  Milstein 
Professor  James  R.  Zeidler 


1992 


TABLE  OF  CONTENTS 


TABLE  OF  CONTENTS . i 

LIST  OF  FIGURES . iv 

LIST  OF  TABLES . ix 

ACKNOWLEDGEMENTS . x 

ABSTRACT . xi 

INTRODUCTION . 1 

1.  MATCHED-FIELD  ARRAY  PROCESSING . 3 

1.1  Introduction . 3 

1.2  Replica  Pressure  Field  Modeling . 3 

1.2.1  Range-Independent  Environment . 4 

1.2. 1.1  Normal  Mode  Model . 5 

1.2. 1.2  Fourier  Integral  Model . 6 

1.2.2  Range-Dependent  Environment . 6 

1.2.2. 1  Normal  Mode  Model . 6 

1. 2.2.2  Parabolic  Equation  Model . 7 

1.3  Beamforming  Processor  Structure . 8 

1.3.1  Bartlett  Method . 8 

1.3.2  Minimum  Variance  Method . 9 

1 .3.3  Eigenvector  Method . 1 1 

1.4  Summary . 12 

2.  1989  SWALLOW  FLOAT  EXPERIMENT . 1 4 

2.1  Introduction . 14 

2.2  Swallow  Float  System  Description . 1 4 

2.3  1989  Experiment  Summary . 16 

2.4  Data  Description . 1 9 

2.4.1  8-KHz  Range  Data . 19 

2.4.2  Acoustic  VLF  Data . 1 9 

2.4.3  Environmental  Data . 22 

2.5  Summary . 25 

3.  SWALLOW  FLOAT  LOCALIZATION . 26 

3.1  Introduction . 26 

3.2  Float  Localization  Techniques . 27 

3.2.1  Least-Squares  Filter . 27 

3.2.2  Kalman  Filter . 27 

3.3  Inputs  to  the  Localization  Filter . 30 

3.3.1  Measurement  Data . 30 

3. 3. 1.1  Travel  Time  Estimates . 30 


i 


3.3. 1.2  Sound  Speed  Estimates . 36 

3.3.2  Measurement  Statistics . 39 

3.3.2. 1  Travel  Time  Variances . 39 

3. 3.2.2  Sound  Speed  Variances . 42 

3.3.3  Initial  Estimates . 42 

3.3.3. 1  Initial  Position  Estimates . 43 

3. 3.3.2  Initial  Velocity  Estimates . 44 

3. 3.3. 3  Initial  Position  and  Velocity  Variance  Estimates . 46 

3.4  Filter  Tuning . 47 

3.4. 1  Least-Squares  Filter . 47 

3.4.2  Kalman  Filter . 48 

3.5  Localization  Results . 50 

3.6  Summary . 56 

4.  SOURCE  LOCALIZATION  WITH  SWALLOW  FLOAT  ARRAY . 57 

4.1  Introduction . 57 

4.2  VLF  Acoustic  Data  Preparation . 58 

4.2.1  Aligning  Float  Time  Bases . 58 

4.2.2  Selecting  Data  Records . 60 

4.2.2. 1  Signal-to-Noise  Ratio . 60 

4.2.2. 2  Spatial  Coherence . 60 

4.2.3  Estimating  the  Array  Covariance  Matrix . 63 

4.2.4  Inverting  the  Array  Covariance  Matrix . 63 

4.3  Acoustic  Propagation  Modeling . 64 

4.3. 1  Environmental  Data . 65 

4.3.2  Modeling  with  Adiabatic  Normal  Modes . 65 

4.3.2. 1  Depth  Eigenfunctions . 67 

4.3.2.2  Horizontal  Wave  Numbers . 70 

4.3.2.3  3-D  Field  Approximation . 71 

4.4  Experimental  Data-Processing  Results . 72 

4.5  Controlled  Simulation . 76 

4.5.1  No  Mismatch  Simulations . 76 

4.5.2  Uncertainty  in  Float  Positions . 80 

4.5.3  Uncertainty  in  Sound  Speed  Structure . 84 

4.6  Summary . 84 

5.  SELF-COHERENT  MATCHED-FIELD  PROCESSING . 89 

5.1  Introduction . 89 

5.2  Literature  Review . 89 

5.2.1  Self-Cohering  Technique . 89 

5.2.2  Focalization  Technique . 91 

5.3  Environment  Adaptation  Technique . 91 

5.3. 1  Global  Optimization  Method . 92 

5.3. 1.1  Motivation . 93 


ii 


5.3. 1.2  Simulated  Annealing  Algorithm . 93 

5.3. 1.3  Fast  Simulated  Annealing  Procedure . 94 

5.3.2  Reducing  the  Parameter  Search  Spaces . 95 

5.4  Processing  Results . 97 

5.4.1  Simulation . 97 

5.4.1. 1  Sound  Speed  Adaptation . 98 

5.4.1.2  Wave  Number  Adaptation . 106 

5.4. 1.3  Remarks . 112 

5.4.2  Experimental  Data . 117 

5.4.2. 1  Sound  Speed  Adaptation . 1 17 

5.4.2.2  Wave  Number  Adaptation . 1 23 

5.5  Summary . 123 

CONCLUSIONS . 128 

APPENDICES . 130 

REFERENCES . 137 


Accession  Fop _ 

NTIS  CRAAI  gT" 

DTIC  TAB  □ 

Unannounced  0 

Ju'-'A  1  f cat.  1  on - - 


By - 

Distribution/ __  _ 

Availability  Codes 
jAvall  and/or 
Dlst  i  [  Special 


LIST  OF  FIGURES 


Figure  1.1  General  beamformer  structure . 9 

Figure  2.1  Schematic  drawing  of  a  typical  Swallow  float . 15 

Figure  2.2  Swallow  floats  deployment  geometry  and  float  retrieval  positions, 

July  1989,  34°  N,  122°  W. . 17 

Figure  2.3  Planned  Swallow  float  deployment  depth.  Also  plotted  in  this 
figure  is  the  sound  speed  profile  based  on  historical  temperature 

and  salinity  data . 1 8 

Figure  2.4  8-kHz  acoustic  range  pulses  received  by  float  1  (freely  drifting) 
as  a  function  of  time  during  the  July,  1989,  experiment.  The 
horizontal  axis  of  the  figure  gives  the  time  in  units  of  Swallow 

float  record  number.  Eighty  records  represent  one  hour. . 20 

Figure  2.5  VLF  acoustic  spectrogram  of  float  1 , 00:28  -  08:28  PST,  9  July  1989. 

The  vertical  axis  of  the  figure  gives  the  time  in  units  of  Swallow 

float  record  number.  Eighty  records  represent  one  hour. . 21 

Figure  2.6  MARK  VI  source  -  Swallow  float  array  geometry,  9  July  1989. 

The  line  of  triangles  marks  the  AXBT  measurements  taken 

between  8  - 10  July  1989 . 23 

Figure  2.7  Sound  speed  profiles  collected  during  the  1989  experiment . 24 

Figure  3.1  Swallow  float  array  geometry  conceptual  diagram . 26 

Figure  3.2  Float  2  listening  to  float  8,  July  1989  experiment.  The  horizontal 
axis  of  the  figure  gives  the  time  in  units  of  Swallow  float  record 

number.  Eighty  records  represent  one  hour. . 32 

Figure  3.3  Float  8  listening  to  float  2,  July  1 989  experiment.  The  horizontal 
axis  of  the  figure  gives  the  time  in  units  of  Swallow  float  record 

number.  Eighty  records  represent  one  hour. . 32 

Figure  3.4  Travel  time  estimates  between  two  freely  drifting  floats  (2  and  8). 

The  horizontal  axis  of  the  figure  gives  the  time  in  units  of  Swallow 

float  record  number.  Eighty  records  represent  one  hour. . 33 

Figure  3.5  Leading  edges  of  each  float’s  detection  of  its  own  surface  echo 
pulses.  The  horizontal  axis  of  the  figure  gives  the  time  in  units  of 

Swallow  float  record  number.  Eighty  records  represent  one  hour. . 34 

Figure  3.6  Float  1 1  surface  echo  pulses,  July  1989  experiment.  The  horizontal 
axis  of  the  figure  gives  the  time  in  units  of  Swallow  float  record 

number.  Eighty  records  represent  one  hour. . 35 

Figure  3.7  Histogram  of  float  1 1  surface  echoes  minus  the  mean  5411  msec 

(binwidth  is  one  msec),  July  1989  experiment . 35 

Figure  3.8  Float  9  listening  to  float  11,  July  1989  experiment.  The  horizontal 
axis  of  the  figure  gives  the  time  in  units  of  Swallow  float  record 
number.  Eighty  records  represent  one  hour. . 37 

iv 


Figure  3.9  Float  1 1  listening  to  float  9,  July  1989  experiment.  The  horizontal 
axis  of  the  figure  gives  the  time  in  units  of  Swallow  float  record 

number.  Eighty  records  represent  one  hour. . 37 

Figure  3.10  Direct  path  travel  time  estimates  between  bottomed  floats  9  and 
11,  calculated  using  constant  mode  values  of  surface  echo  travel 
times,  July  1989  experiment.  The  horizontal  axis  of  the  figure 
gives  the  time  in  units  of  Swallow  float  record  number.  Eighty 

records  represent  one  hour. . 38 

Figure  3. 1 1  Histogram  of  direct  path  travel  time  estimates  minus  the  mean 

4159  msec  (binwidth  is  one  msec),  July  1989  experiment . 38 

Figure  3. 1 2  Travel  time  difference  between  floats  2  and  8  and  a  second-order 
fit,  July  1989  experiment.  The  horizontal  axis  of  the  figure  gives 
the  time  in  units  of  Swallow  float  record  number.  Eighty  records 

represent  one  hour. . 41 

Figure  3.13  Travel  time  difference  (floats  2  and  8)  minus  second-order  fit, 

July  1989  experiment.  The  horizontal  axis  of  the  figure  gives  the 
time  in  units  of  Swallow  float  record  number.  Eighty  records 

represent  one  hour. . 41 

Figure  3. 1 4  Plane  view  of  coordinate  system  used  for  float  localization . 43 

Figure  3.15  Mean  innovation  and  RMS  error  of  the  Kalman  and  least-squares 

filters,  July  1989  experiment . 51 

Figure  3.16  RMS  residual  of  least-squares  filter  position  estimates,  July 
1989  experiment.  The  horizontal  axis  of  the  figure  gives  the 
time  in  units  of  Swallow  float  record  number.  Eighty  records 

represent  one  hour. . 5 1 

Figure  3.17  Float  horizontal  displacement  estimates  using  least-squares 
filter  during  the  July  1989  experiment.  The  circles  mark  the 

starting  positions . 52 

Figure  3.18  Float  depth  estimates  using  least-squares  filter,  July  1989 

experiment.  The  horizontal  axis  of  the  figure  gives  the  time  in 
units  of  Swallow  float  record  number.  Eighty  records  represent 

one  hour. . 53 

Figure  3.19  Float  depth  estimates  using  Kalman  filter,  July  1989  experiment. 

The  horizontal  axes  of  the  figure  give  the  time  in  units  of  Swallow 

float  record  number.  Eighty  records  represent  one  hour. . 54 

Figure  3.20  Float  speed  estimates  derived  from  the  Kalman  filter’s  X,  Y, 

and  Z  velocity  estimates,  July  1989  experiment.  The  horizontal 
axes  of  the  figure  give  the  time  in  units  of  Swallow  float  record 

number.  Eighty  records  represent  one  hour. . 55 

Figure  4.1  Matched-field  processing  block  diagram . 57 

Figure  4,2  (a)  Float  2  listening  to  float  8,  (b)  Float  8  listening  to  float  2, 

(c)  Travel  time  difference  between  floats  2  and  8,  and  (d)  Second- 

order  fit  of  (c) . 59 


v 


Figure  4.3  Acoustic  pressure  power  spectra  for  midwater  floats . 61 

Figure  4.4  (a)  Spatial  coherence  at  14  Hz  and  (b)  Phase  difference  between 

floats  0  and  2  at  14  Hz  during  records  1040  -  1240 . 62 

Figure  4.5  Sound  speed  profiles  between  150°  W  and  122°  W . 66 

Figure  4.6  Sound  speed  profiles  at  the  source  and  Swallow  float  array. . 66 

Figure  4.7  Modal  eigenfunctions  at  frequency  of  14  Hz . . 68 

Figure  4.8  Ray  Daces  for  source  at  90  m . 69 

Figure  4.9  Matched-field  processing  results  using  the  Bartlett  method  for 

the  14-Hz  source  during  record  1145 . 73 

Figure  4.10  Matched-field  processing  results  using  the  MV  method  for  the 

14-Hz  source  during  record  1145 . 74 

Figure  4. 1 1  Matched-field  processing  results  using  the  MUSIC  method  for  the 

14-Hz  source  during  record  1145 . 75 

Figure  4.12  Matched-field  processing  simulation  results  (no  mismatch) 

using  the  Bartlett  method . 77 

Figure  4. 1 3  Matched-field  processing  simulation  results  (no  mismatch) 

using  the  MV  method . 78 

Figure  4.14  Matched-field  processing  simulation  results  (no  mismatch) 

using  the  MUSIC  method . 79 

Figure  4. 15  Matched-field  processing  simulation  of  uncertainty  in  sensor 

positions  using  the  Bartlett  method . 81 

Figure  4. 1 6  Matched-field  processing  simulation  of  uncertainty  in  sensor 

positions  using  the  MV  method . 82 

Figure  4.17  Matched-field  processing  simulation  of  uncertainty  in  sensor 

positions  using  the  MUSIC  method . 83 

Figure  4.18  Matched-field  processing  simulation  of  uncertainty  in  sound 

speed  structure  using  the  Bartlett  method . 85 

Figure  4.19  Matched-field  processing  simulation  of  uncertainty  in  sound 

speed  structure  using  the  MV  method . 86 

Figure  4.20  Matched-field  processing  simulation  of  uncertainty  in  sound 

speed  structure  using  the  MUSIC  method . 87 

Figure  5.1  Environment  adaptation  technique  block  diagram . 92 

Figure  5.2  Matched-field  simulation  (no  mismatch)  of  14-Hz  source  for 
the  three  time  intervals  corresponding  to  the  experimental  data 

collected  during  1989  experiment . 99 

Figure  5.3  Matched-field  simulation  of  environmental  mismatch  for  the 
three  time  intervals  corresponding  to  the  experimental  data 

collected  during  1989  experiment . 100 

Figure  5.4  (a)  Excess  (demeaned)  sound  speed  profiles  computed  from 

30  AXBT  measurements  made  between  140°  W  and  150°  W, 

July  1989,  (b)  Normalized  versions  of  the  first  and  second  EOFs 

derived  from  (a) . 101 


VI 


Figure  5.5  Trajectories  of  the  sound  speed  EOF  coefficients  and  cost  function 

learning  curve  for  a  typical  annealing  run . 103 

Figure  5.6  Energy  surface  as  a  function  of  sound  speed  EOF  coefficients 

computed  by  exhaustive  search  (simulation) . 104 

Figure  5.7  Joint  trajectories  of  sound  speed  EOF  coefficients  for  different 

annealing  runs . . 105 

Figure  5.8  (a)  Sound  speed  profile  derived  from  the  environment  adaptation 

technique  (dotted  line),  “true”  profile  (solid  line),  and  measured 
sound  speed  profile  (dashed  line);  (b)  the  difference  version  of  (a), 
i.e.,  using  the  measured  profile  as  reference  (dashed  line),  the 
“true”  minus  the  measured  and  the  adapted  minus  the  measured 

are  in  solid  and  dotted  lines,  respectively  (simulation) . 107 

Figure  5.9  Self-cohered  source  locations  by  environment  adaptation 

technique  using  sound  speed  EOF  coefficients  as  search  parameters . 108 

Figure  5.10  (a)  Excess  (demeaned)  horizontal  wave  numbers  at  14  Hz 

derived  form  30  AXBT  measurements  made  between  140°  W 
and  150°  W,  July  1989,  (b)  Normalized  versions  of  the  first  and 

second  EOFs  derived  from  (a) . 109 

Figure  5.11  Trajectories  of  the  wave  number  EOF  coefficients  and  cost 

function  learning  curve  for  a  typical  annealing  run  (simulation) . 1 10 

Figure  5.12  Energy  surface  as  a  function  of  wave  number  EOF  coefficients 

computed  by  exhaustive  search  (simulation) . 1 1 1 

Figure  5. 1 3  Joint  trajectories  of  the  wave  number  EOF  coefficients  for 

different  annealing  runs . 1 1 3 

Figure  5.14  Deviations  in  adapted  wave  numbers  and  “true”  wave  numbers 

from  the  measured  model  wave  numbers  (simulation) . 1 14 

Figure  5.15  Self-cohered  source  locations  by  environment  adaptation 
technique  using  wave  number  EOF  coefficients  as  search 

parameters  (simulation) . 115 

Figure  5.16  Range-azimuth  energy  surface  derived  from  repeating  the 

optimization  procedure  for  each  range-azimuth  cell  (simulation) . 116 

Figure  5.17  Matched-field  processing  of  experimental  data  during  the  three 

time  intervals . 118 

Figure  5.18  Range-azimuth  energy  surface  derived  from  repeating  the 

optimization  procedure  for  each  range-azimuth  cell  using  data 

collected  at  01 :47  PST,  July  1989 . 1 1 9 

Figure  5.19  Energy  surface  as  a  function  of  sound  speed  EOF  coefficients 
computed  by  exhaustive  search  corresponds  to  source  location 

during  01:47  PST,  July  1989 . 120 

Figure  5.20  (a)  Sound  speed  profile  (dotted  line)  computed  from  the  optimal 

EOF  coefficients  obtained  from  Figure  5.19  and  the  averaged 
model  sound  speed  profile  (dashed  line),  (b)  The  difference 
version  of  (a)  i.e.,  using  the  measured  profile  as  reference 


vii 


(dashed  line),  the  adapted  minus  the  measured  is  plotted  in 

dotted  line . . . . 121 

Figure  5.2 1  Self-cohered  ‘ounce  locations  by  environmental  adaptation 
technique  3ing  sound  speed  EOF  coefficients  as  the  search 

parameters . 122 

Figure  5.22  Energy  surface  as  a  function  of  wave  number  EOF  coefficients 
computed  by  exhaustive  search  corresponds  to  source  location 

during  01:47  PST,  July  1989 . 124 

Figure  5.23  Adapted  wave  numbers  compared  to  measured  model  wave 

numbers . 125 

Figure  5.24  Self-cohered  source  locations  by  environmental  adaptation 
technique  using  wave  number  EOF  coefficients  as  the  search 
parameters . 126 


LIST  OF  TABLES 


Table  2.1  Planned  and  actual  deployment  depths,  July  1989  experiment . 16 

Table  3. 1  Surface  echo  travel  time  estimates  for  bottomed  floats . 3 1 

Table  3.2  Direct  path  travel  time  estimates  between  bottomed  floats . 36 

Table  3.3  Estimated  variance  of  bottomed  float  travel  time,  July  1989 

experiment . 39 

Table  3.4  Estimated  variance  of  bottomed  float  travel  time  based  on 

predicted  float  movement,  July  1989  experiment . 40 

Table  3.5  Estimated  variance  of  freely  drifting  float  travel  time,  July  1989 

experiment . 42 

Table  3.6  Rough  initial-position  estimate,  record  1003, 00:00  July  9  1989 . 44 

Table  3.7  Initial-position  estimate  produced  by  the  least-squares  filter 

for  record  1003, 00:00,  July  9  1989 . 45 

Table  3.8  Initial  velocity  estimate  (in  meters/second),  record  1003, 00:00, 

July  1989  experiment . 45 

Table  3.9  Estimate  of  initial  position  and  velocity  variance  used  by  the 

Kalman  filter,  July  1989  experiment . 46 

Table  3.10  Bottomed  float  surface  echo  and  direct  path  between  bottomed 

floats  travel  time  estimate,  July  1989  experiment . 48 

Table  3.11  Mean  innovation  power  produced  by  the  Kalman  filter,  July  1989 

experiment . 49 

Table  3.12  Average  float  speed,  July  1989  experiment . 56 

Table  4. 1  Equivalent  source  angles  for  modes  trapped  in  the  water  column 

at  a  frequency  of  14  Hz . 70 

Table  4.2  Eigenvalues  for  modes  trapped  in  the  water  column  at  a 

frequency  of  14  Hz . 71 

Table  4.3  Matched  field  processing  results . 72 

Table  4.4  Matched-field  processing  array  gain . 76 

Table  4.5  Matched-field  simulation  of  sensor  position  errors . 80 

Table  4.6  Matched-field  simulation  of  sound  speed  mismatch . 84 

Table  5. 1  The  five  largest  eigenvalues  obtained  from  eigen-decomposition 

of  the  sound  speed  covariance  matrix . 98 

Table  5.2  Optimization  results  from  nine  runs  using  sound  speed  EOF 

coefficients  as  search  parameters . 02 

Table  5.3  The  five  largest  eigenvalues  obtained  from  eigen -decomposition 

of  the  wave  number  covariance  matrix . ...106 

Table  5.4  Optimization  results  from  nine  runs  using  wave  number  EOF 

coefficients  as  search  parameters . 1 1 2 


ix 


ACKNOWLEDGEMENTS 


I  would  like  to  express  my  sincere  gratitude  to  my  advisor.  Professor  W.  S.  Hodgkiss  for  his 
guidance  and  advice  during  the  research  of  this  thesis.  I  also  wish  to  thank  my  committee  mem¬ 
bers  for  reviewing  this  dissertation  and  offering  constructive  criticism. 

I  am  indebted  to  the  graduate  academic  fellowship  program  sponsored  by  the  Naval  Com¬ 
mand,  Control  and  Ocean  Surveillance  Center’s  RDT&E  Division  (formally  Naval  Ocean  Sys¬ 
tems  Center)  and  the  support  of  the  N03A  Block  program  of  the  Office  of  Naval  Technology. 
Without  their  support,  my  graduate  studies  at  UCSD  would  not  have  been  possible.  Special  thanks 
are  due  to  Dr.  C.  E.  Persons  Dr.  D.  K.  Barbour,  and  Dr.  J.  A.  Roese  of  NRaD  for  their  continuous 
encouragement  and  unwavering  support  over  the  last  few  years. 

During  a  major  portion  of  this  research,  the  Marine  Physical  Laboratory  (MPL)  kindly  provid¬ 
ed  the  office  space  and  computing  facilities.  The  hospitality  of  MPL  during  this  period  is  greatly 
appreciated.  I  have  enjoyed  sharing  ideas  and  learning  from  several  members  of  Professor  W.  S. 
Hodgkiss’  acoustic  signal  processing  research  group  at  MPL.  Special  thanks  go  to  Dr.  G.  L. 
D’Spain  for  sharing  his  expertise  with  the  Swallow  iloats  and  Dr.  J.  M.  Tran,  a  former  member  of 
MPL,  for  many  insightful  discussions  on  array  processing.  I  also  thank  the  members  of  the  Swal¬ 
low  float  team  for  preparing  the  sea  experiment  and  collecting  the  experimental  data  needed  in 
this  dissertation 

Finally,  I  am  grateful  to  my  family  for  their  support  and  for  enduring  the  countless  hours  of 
my  graduate  studies  over  these  last  few  years.  It  is  to  them  I  dedicate  this  dissertation. 

The  Swallow  Float  experiment  was  sponsored  by  the  Office  of  Naval 
Technology  under  Naval  Research  Laboratory  Contract  N00014-88-K-2040. 


x 


ABSTRACT  OF  THE  DISSERTATION 


VLF  Source  Localization  with  a  Freely  Drifting  Sensor  Array 

by 


George  Chia-Jen  Chen 

Doctor  of  Philosophy  in  Electrical  Engineering 
(Applied  Ocean  Science) 

University  of  California,  San  Diego,  1992 
Professor  William  S.  Hodgkiss,  Jr.,  Chair 


Source  localization  and  tracking  capability  of  the  freely  drifting  Swallow  float  volumetric  ar¬ 
ray  is  demonstrated  with  matched-field  processing  (MFP)  technique  using  data  collected  during 
the  1989  Swallow  float  experiment  conducted  in  the  Northeast  Pacific.  Marine  Physical  Laborato¬ 
ry’s  set  of  nine  freely  drifting,  infrasonic  sensors,  capable  of  recording  ocean  ambient  noise  in  the 
1-  to  25-Hz  range,  was  deployed  to  span  the  water  column  of  4100-m  depth,  with  horizontal  aper¬ 
ture  on  the  order  of  5  km.  Even  though  the  floats  were  freely  drifting,  their  positions  were  deter¬ 
mined  with  the  8-kHz  internal  navigation  system  and  a  postprocessing  least-squares-based 
localization  procedure.  The  rms  position  error  estimated  by  the  float  localization  procedure  is  less 
than  5  meters,  which  is  within  the  desired  accuracy  of  one-tenth  of  a  wavelength  at  the  highest 
frequency  of  interest,  25  Hz  (6  m),  in  order  to  effectively  beamform  the  acoustic  data.  Analysis  of 
the  acoustic  data  showed  high  signal-to-noise  ratio  and  high  magnitude-squared  coherence  at  14 
Hz  among  all  floats  during  some  time  intervals.  The  1 4-Hz  was  a  continuous  wave  (cw)  tone  pro¬ 
jected  by  a  source  of  opportunity  deployed  about  2500  km  from  the  Swallow  float  array.  The  high 
coherence  among  the  floats  provided  an  opportunity  to  matched-field  process  the  acoustic  data. 
The  replica  pressure  field  was  modeled  with  an  adiabatic  normal-mode  numerical  technique  using 
the  environ  mental  data  collected  during  the  experiment.  Initial  MFP  of  the  experimental  data  re¬ 
vealed  difficulties  in  estimating  the  source  depth  and  range  while  the  source  azimuth  estimate  was 
quite  successful.  The  main  cause  of  the  MFP  performance  degradation  was  incomplete  knowl¬ 
edge  of  the  environment.  An  environment  adaptation  technique  using  a  global  optimization  algo¬ 
rithm  was  developed  to  alleviate  the  environmental  mismatch  problem.  With  limited  knowledge 
of  the  environment  and  a  known  location  of  the  14-Hz  source  during  a  selected  time  interval  ac¬ 
cording  to  the  source  log,  the  ocean-acoustic  environment  can  be  adapted  to  the  acoustic  data  in  a 
matched-field  sense.  Using  the  adapted  environment,  the  14-Hz  source  was  successfully  localized 
and  tracked  in  azimuth  and  range  within  a  region  of  interest  using  the  MPT  technique  for  a  period 
of  6  hours. 
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INTRODUCTION 


Traditionally,  source  localization  has  relied  on  the  processing  of  assumed  plane-wave  fronts 
received  by  spatially  distributed  sensors  to  estimate  the  source  bearing  or  vertical  angle  of  arriv¬ 
als.  In  reality,  the  ocean  acoustic  channel  is  extremely  complex  due  to  refractive  and  multipath  ef¬ 
fects.  Assumption  of  plane-wave  arrivals  in  the  processing  scheme  can  lead  to  severe  degradation 
of  the  estimate.  Matched-field  processing  (MFP)  has  been  proposed  [Bucker,  1976]  to  actually 
use  the  complex  ocean  acoustic  properties  to  improve  source  detection  and  localization.  MFP  in¬ 
volves  the  correlation  of  the  actual  acoustic  pressure  field  measured  at  the  array  with  a  predicted 
field  due  to  a  source  at  an  assumed  location  deriving  from  an  acoustic  propagation  model.  A  high 
degree  of  correlation  between  the  measured  field  and  the  predicted  field  indicates  a  likely  source 
location.  Matched-field  processing  of  the  acoustic  wavefield  has  shown  that  when  sufficient  envi¬ 
ronmental  characterizations  (e.g.,  sound  speed  profile,  bathymetry,  sediment  properties)  are  avail¬ 
able,  rather  remarkable  detection  and  estimation  results  can  be  obtained.  Most  available  matched- 
field  work  has  been  for  rather  simple  propagation  situations  (e.g.,  range-independent  environment 
or  shallow-water  environment)  and  much  of  the  work  has  been  restricted  to  vertical-line  arrays 
[Bucker,  1976;  Porter  et.  al„  1987;  Fizell,  1987;  Baggeroeret.  al„  1988;  Ozard,  1989;  Feuillade 
et.  al.,  1990;  Zala  and  Ozard,  1990;  Sotirin  et.  al.,  1990;  Tran  and  Hodgkiss,  1991]. 

Although  matched-field  processing  offers  an  appealing  approach  to  the  underwater  source  de¬ 
tection  and  estimation  problem,  a  common  difficulty  with  this  technique  occurs  when  the  environ¬ 
ment  information  is  inaccurate.  A  “mismatch”  occurs  between  the  measured  data  and  the 
modelled  pressure  field,  and  the  performance  of  the  MFP  is  degraded  and  leads  to  errors  in  the  es¬ 
timation  of  the  source  location  [Tolstoy,  1989;  Feuillade  et.  al.,  1989;  Hamson  and  Heitmeyer, 
1989;  Gingras,  1989;  Daugherty  and  Lynch,  1990], 

The  focus  of  this  dissertation  research  is  twofold:  (1)  to  demonstrate  the  match-field  source  lo¬ 
calization  and  tracking  capability  of  the  Swallow  float  freely  drifting  volumetric  array  using  ex¬ 
perimental  data,  and  (2)  to  study  the  environmental  mismatch  problem  and  investigate  the 
technique  that  may  neutralize  the  effect  caused  by  imprecise  knowledge  of  the  environment  and 
thereby  lead  to  MFP  localization  performance  enhancement. 

This  dissertation  is  organized  as  follows.  Chapter  1  presents  an  overview  of  the  most  com¬ 
monly  used  low-frequency  acoustic  propagation  models  and  array  processing  power  estimators. 
Five  acoustic  propagation  models  are  presented  in  categories  reflecting  the  ability  of  the  model  to 
handle  environmental  range  dependence.  Three  array  processing  methods  (Bartlett,  Minimum 
Variance,  and  MUSIC)  are  reviewed  in  the  context  of  matched-field  processing,  and  the  basic 
concepts  and  properties  associated  with  each  method  are  discussed. 
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Chapter  2  gives  a  brief  description  of  the  Swallow  float  system  and  a  summary  of  the  July 
1989  Swallow  float  experiment.  Data  collected  during  the  experiment,  which  includes  8-kHz 
range  data,  VLF  acoustic  data,  and  environmental  data,  are  presented. 

In  chapter  3,  the  8-kHz  range  data  collected  by  the  Swallow  floats  during  the  1989  experi¬ 
ment  are  used  to  estimate  float  positions  as  a  function  of  time.  Two  methods  (least-squares  and 
Kalman  filters)  used  in  localizing  the  floats  are  reviewed,  and  their  results  are  compared. 

Chapter  4  presents  the  results  of  matched-field  processing  on  the  14-Hz  continuous  wave  (cw) 
tone  collected  by  the  Swallow  floats  during  the  1989  experiment.  Issues  related  to  float  time-base 
alignment,  data  selection  criteria,  and  array  covariance  matrix  estimation  are  addressed.  An  adia¬ 
batic  normal-mode  modeling  technique  and  the  environmental  data  collected  during  the  experi¬ 
ment  are  used  to  model  the  acoustic  replica  pressure  fields.  All  three  array  processors  are 
implemented  and  used  to  compute  the  matched-field,  range-depth  and  range-azimuth  ambiguity 
surfaces  on  the  experiment  data.  Controlled  simulation  is  also  performed  to  aid  in  interpreting  the 
experimental  data  processing  result. 

Chapter  5  addresses  the  environmental  mismatch  problem.  Two  techniques,  self-cohering  and 
focalization,  reported  in  the  literature  are  briefly  reviewed.  An  environment  adaptation  technique 
using  a  global  optimization  procedure  is  proposed  to  enhance  the  MFP  localization  performance. 
The  optimization  procedure  is  implemented  using  the  empirical  orthogonal  function  (EOF)  ap¬ 
proach  to  reduce  the  parameter  spaces  and  a  fast  simulated  annealing  algorithm  to  search  for  the 
optimal  environmental  parameter  values.  Two  types  of  environmental  parameters  are  considered: 
sound  speed  structure  and  modal  wave  number.  The  technique  is  illustrated  using  both  simulation 
and  experimental  data. 

Lastly,  a  summary  of  this  dissertation  and  a  discussion  of  area  of  future  research  is  given. 
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1.  MATCHED-FIELD  ARRAY  PROCESSING 


1.1  Introduction 

The  problem  of  coherently  summing  the  outputs  from  a  number  of  spatially  distributed  sen¬ 
sors  to  form  a  spatial  filter  is  known  as  beamforming  [De  Fatta,  Lucas,  and  Hodgkiss,  1988],  In 
underwater  signal  processing,  conventional  beamforming  can  be  viewed  as  correlating  the  pres¬ 
sure  field  measured  at  the  sensors  with  a  plane-wave  replica  at  an  assumed  incidence  angle  on  the 
array.  A  peak  in  the  plane  of  azimuthal  and  elevation  angle  indicates  the  presence  and  direction  of 
the  source.  However,  due  to  the  refractive  and  multipath  effects  in  the  ocean,  the  acoustic  pressure 
field  received  by  an  array  from  a  low-frequency  narrowband  point  source  cannot  be  expected  to 
be  planar.  The  nature  of  the  field  is  a  function  of  range,  depth,  and  azimuth  of  the  source  as  well  as 
a  function  of  the  sound  speed  structure  and  the  characteristics  of  the  ocean  surface  and  bottom.  If 
array  detection  and  localization  performance  are  to  be  maximized,  the  assumed  replica  for  the  sig¬ 
nal  field  must  match  the  actual  field  as  closely  as  possible.  Then  the  steering  vector  must  scan  in 
range,  depth,  and  azimuth,  rather  than  just  in  directions.  Scanning  in  range,  depth,  and  azimuth  re¬ 
quires  a  prediction  of  the  spatial  structure  of  the  received  acoustic  field  for  each  potential  source 
location.  The  resulting  spatial-filtering  process  is  referred  to  as  matched-field  processing. 
Matched-field  processing  is  a  generalization  of  conventional  plane-wave  beamforming.  In  es¬ 
sence,  MFP  matches  the  measured  field  at  the  array  with  a  realistic  acoustic  model  prediction  of 
the  pressure  due  to  a  source  at  an  assumed  location.  In  this  case,  a  peak  in  the  range/depth/azi¬ 
muth  ambiguity  surfaces  indicates  the  presence  and  location  of  the  source. 

To  allow  source  detection  and  localization  in  the  matched-field  sense,  the  correlation  between 
the  measured  pressure  field  and  the  replica  field  must  be  measured.  Thus,  the  major  signal  pro¬ 
cessing  functions  for  matched-field  processing  include  acoustic  propagation  modeling,  which  pre¬ 
dicts  the  replica  pressure  field,  and  beamformer  power  estimation,  which  represents  the 
agreement  between  actual  and  replica  fields. 


1.2  Replica  Pressure  Field  Modeling 

Formulations  of  acoustic  propagation  modeling  generally  begin  with  the  wave  equation.  The 
wave  equation  is  a  partial  differential  equation  relating  the  acoustic  pressure  or  the  velocity  poten¬ 
tial  to  the  space  and  time,  and  may  be  written  as 


V2d>  -  — — 1 - 

c  ( x,y,z ) 


d? 


(1.1) 


where  V2  is  the  Laplacian  operator,  d>  is  the  velocity  potential  function,  and  c  is  the  speed  of 
sound.  If  harmonic  solution  is  assumed  for  the  potential  function 
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(1.2) 


O.  -iwf 

=  4>  e 


where  4>  is  the  time-independent  potential  function,  and  o)  is  the  source  frequency. 
Then  Equation  (1.1)  can  be  reduced  to 


V2<J>  -I-  A:24>  =  0 


(1.3) 


or  in  cylindrical  coordinates, 


O) 

where  k  =  — . 

c 

Equation  (1 .3)  is  the  time-independent  (or  frequency  domain  wave  equation)  also  known  as  the 
Helmholtz  equation.  Various  theoretical  approaches  are  applicable  to  the  Helmholtz  equation  de¬ 
pending  upon  the  specific  geometrical  assumptions  made  for  the  propagation  and  the  type  of  solu¬ 
tion  chosen  for  <}>.  In  practice,  the  Helmholtz  equation  is  too  complex  for  analytical  solutions; 
numerical  methods  are  generally  used. 

To  describe  the  different  wave  theory  numerical  modeling  approaches,  it  is  convenient  to  di¬ 
vide  them  according  to  range-independent  and  range-dependent  types.  Range  independence 
means  that  the  model  assumes  a  cylindrical  symmetry  for  the  environment  (i.e.,  a  horizontally 
stratified  ocean  in  which  properties  vary  only  as  a  function  of  depth).  Range  dependence  indicates 
that  some  properties  of  the  ocean  medium  are  allowed  to  vary  as  a  function  of  range  (r)  and  azi¬ 
muth  (0)  between  the  source  and  sensor,  in  addition  to  a  depth  ( z )  dependence.  Such  range-vary¬ 
ing  properties  commonly  include  sound  speed  structure,  bathymetry,  and  sediment  properties.  The 
following  sections  review  low-frequency  acoustic  modeling  techniques  commonly  referenced  in 
the  literature.  The  presentation  closely  follows  an  excellent  paper  [Jensen,  1988]  and  a  recent  text 
[Etter,  1991]. 

1.2.1  Range-Independent  Environment 

If  the  ocean  is  assumed  to  be  perfectly  horizontally  stratified  and  azimuthally  symmetrical, 
i.e.,  the  sound  speed  c  (z)  varies  only  with  depth,  the  solution  for  the  potential  function  <{>  can  be 
written  as  the  product  of  a  depth  function  Z(z)  and  a  range  function  R(r)  using  the  separation  of 
variables  technique: 


92<l>  1 3<1>  3  2<b  , 


(1.4) 


<D  =  Z (z)R  (r) 


(1.5) 
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If  the  separation  constant  is  designated  £2,  the  separated  equations  are  [Boyles,  1984] 

?Jt+(k2-¥)Z  =  0  (1.6) 

dz2 

^  +  -f?  +  |2tf  =  0  (1.7) 

dr 2  rdr 

Equation  (1.6)  is  the  depth  equation,  known  as  normal  mode  equation,  which  describes  the  stand¬ 
ing  wave  portion  of  the  solution;  its  solution  is  the  Green’s  function.  Equation  (1.7)  is  the  range 
equation,  a  zero-order  Bessel  equation,  which  describes  the  traveling  wave  portion  of  the  solu¬ 
tion;  its  solution  can  be  written  in  terms  of  a  zero-order  Hankel  function.  The  full  solution  for  <{> 
can  then  be  expressed  by  an  infinite  integral,  assuming  a  single-frequency  point  source: 


4>{r,z)  =  \G{z,z0\l)H^  (gr)lc^  (1.8) 

■— OO 

where  G( )  is  the  depth-dependent  Green’s  function,  and  Hq1\  )  is  the  zero-order  Hankel  func¬ 
tion  for  the  first  kind,  and  z0  is  the  source  depth. 

1.2.1.1  Normal  Mode  Model 

To  obtain  the  normal  mode  solution  to  the  wave  equation,  the  Green’s  function  is  expanded 
into  a  complete  set  of  normal  mode  functions  um.  The  eigenvalues  (or  the  horizontal  wave  num¬ 
bers),  which  are  the  resulting  value  of  the  separation  constants,  are  represented  by  'E>m.  These 
eigenvalues  represent  the  discrete  set  of  values  for  which  solutions  of  the  modal  eigenfunctions 
um  exist.  The  spectrum  of  eigenvalues  generally  consists  of  both  a  discrete  part  and  a  continuous 
part.  By  ignoring  the  continuous  part,  which  correspond  to  sound  interacting  with  the  bottom  at 
angles  above  the  critical  grazing  angle  (or  nonpropagating  source  nearfield)  and  by  replacing  the 
Hankel  function  expression  by  its  far-field  asymptotic  expansion  for  \r  »  1 ,  the  solution  for  4>  can 
be  written  as  a  modal  sum: 


M 

<$>(r,z)  =  const  £  u m(z0)um(z) 

m  =  1 


(1.9) 


where  M  is  the  number  of  propagating  modes. 


1.2.1.2  Fourier  Integral  Model 

An  alternative  approach  to  (1.8)  is  to  replace  the  Hankel  function  expression  in  (1.8)  by  the 
first  term  in  the  asymptotic  expansion: 


<Maz) 


(1.10) 


The  infinite  integral  is  then  evaluated  by  means  of  the  fast  Fourier  transform,  which  provides  val¬ 
ue  of  the  potential  function  <J>  at  n  discrete  range  points  for  a  given  source-receiver  geometry.  This 
approach  requires  greater  computation  effort  than  the  normal  mode  method  but  does  include  the 
nonpropagating  source  nearfield. 


1,2.2  Range-Dependent  Environment 


1. 2.2.1  Normal  Mode  Model 

As  discussed  in  the  last  paragraph,  the  normal  mode  model  assumes  range  independence.  Ex¬ 
tension  to  range  dependence  can  be  accomplished  either  by  mode  coupling,  which  considers  the 
energy  scattered  from  a  given  mode  into  other  modes,  or  by  invoicing  the  adiabatic  assumption, 
which  assumes  that  all  energy  in  a  given  mode  transfers  to  the  corresponding  mode  in  the  new  en¬ 
vironment,  provided  that  the  environmental  variations  in  range  are  gradual.  This  section  considers 
two  such  extensions. 

Coupled  Mode  Method 

For  range-dependent  acoustic  propagation  environment,  the  range  between  source  and  sensor 
can  be  divided  into  a  number  of  range  regions,  each  with  range-invariant  properties  but  with  al¬ 
lowance  for  arbitrary  variation  of  environmental  parameters  with  depth  within  each  region.  Since 
the  range  dependence  takes  place  discretely  at  the  interfaces  between  regions,  the  solution  within 
a  range-independent  region  can  be  constructed  using  the  standard  normal  mode  solution  and  inter¬ 
face  conditions.  A  complete  two-way  solution  for  wave  propagation  can  be  written  as  a  sum  of  lo¬ 
cal  modes  with  unknown  excitation  coefficients  determined  by  requiring  continuity  of  pressure 
and  horizontal  particle  velocity  across  region  boundaries  [Evans,  1983],  Computationally,  this  ap¬ 
proach  requires  the  solution  of  a  whole  family  of  normal  mode  problems,  one  for  each  region,  fol¬ 
lowed  by  a  large  banded,  block  linear  system.  The  coupled  mode  technique  does  require 
considerable  computational  power  but  provides  complete  wave  solutions  with  backscattering  in¬ 
cluded. 
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Adiabatic  Approximation  Method 


If  the  cross-mode  coupling  in  the  coupled  mode  method  is  neglected,  we  can  derive  a  simple 
first-order  approximation  to  the  full-coupling  problem,  called  the  adiabatic  normal  mode  solution. 
It  is  assumed  that  the  energy  c  stained  in  a  given  mode  at  the  source  range  will  stay  in  that  partic¬ 
ular  mode  throughout  the  range-varying  environment  and  not  couple  to  neighboring  modes.  This 
assumption  is  accurate  for  a  slow-varying,  range-dependent  medium.  The  solution  for  0  in  a  weak 
range-dependent  environment  can  be  expressed  as  [Pierce,  1965;  Brekhovskikh  and  Lysanov, 
1982] 


M 

0(r,z)  =  const  £  um(z0;0)um(z;r) 

m  =  1 


(i.ii) 


The  adiabatic  mode  technique  is  computational  efficient  since  no  mode  coupling  effects  need  to 
be  considered. 


1.2.2.2  Par^olic  Equation  Model 

The  parabolic  equation  (PE)  approach  is  a  far-field,  narrow-angle  approximation  to  the  wave 
equation  [Tapped,  1977].  An  approximate  solution  to  2-D  propagation  problems  can  be  obtained 
from  the  parabolic  equation,  which  is  derived  by  assuming  a  field  solution  of  the  form: 


0(r,z)  =  v(r,z) 


exp  ( ik0r ) 

~F~ 


(1.12) 


where  k0  is  an  average  horizontal  wavenumber  of  the  propagating  wave  spectrum.  By  substitut¬ 
ing  the  expression  for  0  into  (1.1),  the  equation  for  \j/  (r,  z)  can  be  simplified  to 


afy  aV 
a7  +  a? 


+  +*o 


(n2-  1)V  =  o 


(1.13) 


where  n  =  k(r,z)/  ka  is  the  refraction  index.  Further  assume 


aV 

dr 


llr  ^ 

2  ‘  2k°dr 


0.14) 


which  is  the  paraxial  approximation  and  is  equivalent  to  saying  that  if  the  acoustic  field  were  rep¬ 
resented  by  rays,  these  would  be  inclined  only  at  small  angles  to  the  horizontal.  The  paraxial  ap- 
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proximation  is  the  heart  of  the  parabolic  equation  method,  and  by  combining  (1.13)  and  (1.14)  we 
get 


?y+2ikl¥  +  kl(n2-\)q>  =  0  (1.15) 

dz  dr 

This  is  the  standard  parabolic  equation  that  represents  a  wave  propagating  outwards  and  close  to 
the  horizontal.  This  approximation  takes  us  from  a  boundary-value,  wave-equation  problem  to  an 
initial-value,  parabolic  equation  that  leads  itself  to  a  marching-solution  technique  when  the  initial 
condition  is  known.  Note  that  PE  is  an  approximation  technique  that  neglects  backscattering  and, 

o 

more  importantly,  assumes  energy  to  propagate  within  a  limited  angular  spectrum  of  ±40  with  re¬ 
spect  to  the  horizontal. 

1.3  Beamforming  Processor  Structure 

Conventional  beamforming  procedures  use  plane-wave  steering  vectors  as  replicas  of  the  true 
acoustic  pressure  field  and  estimate  the  source  bearing  by  finding  the  maximum  beamformer  out¬ 
put  that  represents  the  agreement  between  actual  and  replica  fields.  The  matched-field  processor 
structure  takes  the  same  form  as  the  conventional  processor  except  that  the  plane-wave  replica 
vector  is  replaced  by  a  replica  vector  derived  from  an  acoustic  propagation  model  for  the  oceanic 
waveguide  of  interest.  Three  such  power  estimators  commonly  used  in  conventional  plane-wave 
array  processing  that  differ  in  their  resolution  are  reviewed  in  the  context  of  matched-field  pro¬ 
cessing. 

1.3.1  Bartlett  Method 

Consider  the  general  beamformer  structure  shown  in  Figure  1.1 .  The  array  of  interest  consists 
of  M  sensors  of  known  but  arbitrary  geometry.  With  xi  (/)  representing  the  output  at  the  ich  sen¬ 
sor  in  frequency  domain  and  wi  representing  the  corresponding  desired  frequency  domain 
weighting  factor,  the  beamformer  output  can  be  written  as  matrix  form: 

y(/)=WwX(/)  (1.16) 

The  corresponding  beamformer  output  power  is  given  by  [Pillai,  1989] 

P  =  E[m2]  =  E  [W^XX^W]  =  WwRW  (1.17) 

where  R  =  E  [XXW]  is  the  cross  spectral  density  matrix  or  array  covariance  of  the  sensor  out¬ 
puts. 


8 


If  we  replace  the  weight  vector  W  by  E  (r,  z,  0) ,  the  replica  pressure  field  vector  due  to  a  narrow- 
band  source  at  (r,  z,  0) ,  then  the  output  of  the  beamformer  can  be  expressed  as 

P Barden  =  EH  (r,  z,  0)  R  E  (r,  z,  0)  (1.18) 

Since  (1.18)  has  the  same  form  as  the  Bartlett  estimate  in  spectral  estimation,  (1.18)  is  also  called 
a  Bartlett  processor  in  array  processing  literature.  The  Bartlett  processor  is  known  to  be  robust  but 
offers  low  resolution. 


Figure  1.1  General  beamformer  structure. 


1.3.2  Minimum  Variance  Method 

Recently,  studies  of  matched-field  processing  have  also  centered  upon  the  use  of  the  maxi¬ 
mum-likelihood  method  by  Capon  [1969]  (also  known  as  the  minimum  variance  [MV]  method) 
for  its  sidelobe  suppression  capability  and  high  resolution  of  the  mainlobe.  Basically,  the  algo¬ 
rithm  adaptively  constructs  a  weight  vector  W  that  yields  a  minimum  mean-square  output 
WwRW  to  the  noise  field,  with  the  constraint  that  signal  field  is  passed  with  unity  gain  WwE  =  1 
when  the  postulated  source  location  is  equal  to  the  true  source  location: 
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(1.19) 


minimize  WwR  W 

VV 

subject  to  WhE  =  1  (1.20) 

In  other  words,  the  weight  vector  minimizes  the  quadratic  function  [Kanasewich,  1975]: 


W^RW+MW^E-l]  (1.21) 

w 

where  X  is  an  undetermined  Lagrange  multiplier. 

If  we  express  (1.21)  as 

^{W"R  W  +  X[WwE- 1]>  =  0  (1.22) 

(1.22)  can  be  easily  reduced  to 

RW  =  ^E  (1.23) 

Multiplying  by  the  inverse  of  R  yields 

W  =  ^R_,E  (1.24) 

To  determine  the  Lagrange  multiplier,  X,  substitute  W  into  the  constraint  equation  (1.20). 


X  = 


-2 


EWR-1E 


The  optimal  weight  vector  is  derived  as 


W  = 


R~*E 


EWR_1E 

Substituting  (1.26)  into  (1.17),  the  corresponding  beamformer  output  power  is 

1 


P MV  (r»  B)  — 


E//R1E 


(1.25) 


(1.26) 


(1.27) 
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In  general,  the  Bartlett  beamformer  has  significant  sidelobes  that  often  are  indistinguishable 
from  the  mainlobe.  The  MV  matched-field  processor  has  better  suppression  of  the  sidelobes  and 
high  resolution  of  the  mainlobe,  but  it  is  more  sensitive  to  errors  in  replica  vectors. 

1.3.3  Eigenvector  Method 

Intuitively,  it  appears  that  we  could  further  improve  the  resolution  properties  of  the  minimum 
variance  beamformer  if  the  denominator  of  (1.27)  could  be  made  smaller  in  those  locations  that 
correspond  to  signals  and  larger  in  locations  corresponding  to  noise.  The  eigenvector  beamformer 
is  derived  based  on  this  intuition  [Johnson,  1982].  Assuming  we  have  M  sensors  and  K  sources,  it 
is  well  known  that  the  K  largest  eigenvalues  of  the  array  covariance  matrix  R  are  associated  with 
the  sources.  The  corresponding  K  eigenvectors  span  the  signal  subspace,  that  is,  they  are  linear 
combinations  of  the  signal  vectors.  Furthermore,  the  M-K  remaining  eigenvectors  span  a  noise 
subspace  orthogonal  to  the  signal  subspace.  We  can  express  R  in  terms  of  its  eigenvalues  X  and 
eigenvectors  V(,  i  -  1, ...,  M,  using  the  spectral  theorem  fora  matrix: 

M 

R=XW'f  o-28> 

i  =  I 

The  inverse  of  R  has  an  equally  simple  form: 


« =  l  ' 


Next,  we  form  a  truncated  or  noise  covariance  matrix  by  including  only  the  M-K  eigenvectors  that 
span  the  noise  subspace.  This  matrix  is  denoted  by  R^  and  the  inverse  of  RN  can  be  expressed  as 

-  I  rv/vf  0-30) 

i  =  1  < 

The  noise  covariance  then  is  used  in  (1.26)  to  find  the  corresponding  optimum  weight  vector  for 
the  eigenvector  beamformer: 


W  = 


R*E 

E"R"lE 


The  corresponding  eigenvector  beamformer  output  power  is 


(1.31) 
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1 


(1.32) 


P EV  (A  Z,  0)  — 


E^R'1  E 


M-K 


Ir 

i  =  i  ' 


e7/v.  2 


Because  vectors  in  the  signal  space  are  orthogonal  to  the  noise  subspace  spanned  by  V/t 
i  =  1, M  -  K,  the  inner  product  EWV(.  is  zero  if  E  is  in  the  signal  subspace.  Thus  with  the  as¬ 
sumption  of  uncorrelated  noise,  the  denominator  of  (1.32)  is  equal  to  zero  for  source  locations. 
The  power  at  the  output  of  the  beamformer  given  by  (1.32)  therefore  is  infinite  in  those  locations 
in  the  signal  subspace. 


A  related  eigenvector  method  called  MUSIC  (Multiple  Signal  Classification)  [Schmidt,  1 986] 
sets  the  small  eigenvalues  equal  to  unity  in  the  formation  of  the  noise  covariance  matrix: 


M-K 

(MUSIC)  ~  X  vy?  (1.33) 

i  =  1 

Equation  (1.32)  can  be  rewritten  as 


P MUSIC  (r’Z’Q)  -  w/rTi  ~  -  M^K^ -  (L34) 


The  effect  of  setting  the  small  eigenvalues  to  unity  is  to  “whitten”  those  locations  that  do  not  cor¬ 
respond  to  the  source  signals.  The  resolution  capability  of  the  MUSIC  method  is  similar  to  that  of 
the  eigenvector  method. 

1.4  Summary 

Matched-field  processing  is  a  generalization  of  plane-wave  beamforming  in  which  the  plane- 
wave  replicas  are  replaced  by  solutions  to  the  wave  equation  for  the  environment  of  interest. 
There  are  two  stages  in  applying  MFP  to  narrowband  cw  data  received  by  the  array:  full-wave 
field  modeling  and  power  estimation.  Full-wave  field  modeling  involves  the  use  of  an  acoustic 
propagation  numerical  model  for  predicting  the  acoustic  pressure  field  in  a  waveguide  of  known 
environmental  properties.  Based  on  the  model,  the  acoustic  replica  pressure  field  at  each  sensor  of 
the  array  may  be  computed  for  a  source  at  a  particular  position  with  respect  to  the  array.  Power  es¬ 
timation  consists  of  matching  the  replica  pressure  fields  with  the  measured  data.  Peak  in  the  pow¬ 
er  output  with  respect  to  range,  azimuth,  and  depth  may  be  taken  as  estimates  of  the  position  of 
the  source  within  the  waveguide.  In  this  chapter,  five  common  low-frequency  acoustic  propaga¬ 
tion  modeling  techniques  were  reviewed.  Area  of  applicability  for  each  modeling  technique  was 
described.  Three  of  the  most  commonly  used  array  processing  techniques  were  reviewed  in  the 
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2.  1989  SWALLOW  FLOAT  EXPERIMENT 

2.1  Introduction 

Data  presented  in  this  thesis  were  collected  by  the  Marine  Physical  Laboratory’s  (MPL’s) 
Swallow  floats  between  8  and  9  July  1989,  approximately  150  km  w-est-northwest  of  Pt.  Arguello, 
California,  near  34°  50'  N,  122°  20'  W,  A  detailed  trip  report  for  the  experiment  is  contained  in 
[Chen  et.  al.,  1990].  This  chapter  describes  the  Swallow  float  system  description;  a  summary-  of 
the  experiment;  the  Swallow  float  system  description;  and  representative  data  collect  i  during  the 
experiment. 

2.2  Swallow  Float  System  Description 

Over  the  last  several  years,  MPL  has  designed  and  developed  a  set  of  12  acoustic  sensors  that 
are  neutrally  buoyant  and  freely  drifting  when  deployed  in  the  ocean.  The  sensors  are  called  Swal¬ 
low  floats  in  honor  of  J.  C.  Swallow  who  used  freely  drifting  floats  to  measure  deep  ocean  cur¬ 
rents.  The  MPL  Swallow  floats  are  used  to  measure  acoustic  energy  in  the  very  low  frequency 
(VLF)  band,  which  is  generally  defined  as  1  to  25  Hz.  The  Swallow  floats  are  designed  to  operate 
without  tether  so  that  their  measurements  are  not  contaminated  by  tether  strumming  noise  and 
flow  noise. 

As  illustrated  in  Figure  2. 1 ,  the  Swallow  float  system  is  a  1 7-inch-diameter  glass  float  contain¬ 
ing  three  orthogonally  oriented  geophones  used  as  the  acoustic  particle  velocity  sensor,  a  com¬ 
pass,  the  electronics  and  power  supply.  External  to  the  sphere  are  a  hydrophone  for  measuring 
VLF  acoustic  pressure,  an  8-kHz  acoustic  transducer  for  transmitting  and  receiving  acoustic  rang¬ 
ing  signals,  an  optical  flash,  a  radio  beacon,  and  the  release  ballast  [D’Spain  et.  al.,  1991]. 

In  operation,  the  floats  are  deployed  and  ballasted  to  neutral  buoyancy  at  the  desired  depths. 
While  deployed,  each  float  records  signals  from  the  three  geophones  and  the  hydrophone  sampled 
at  50  Hz,  the  compass,  and  the  8-kHz  range  pulse  arrival  times.  Each  float  continuously  fills  a  12- 
kB  buffer  memory,  then  periodically  writes  this  buffer  to  tape.  Float  time  is  divided  into  45-sec- 
ond  periods  consisting  of  44  seconds  during  which  the  geophone  and  hydrophone  outputs  are 
sampled  and  stored  in  the  buffer  and  1  second  during  which  the  buffer  contents  are  written  to  tape. 
No  data  are  recorded  during  the  1 -second  data  dump.  Limited  by  tape-recorder  capacity  (i.e.,  17 
Mbytes)  the  data  acquisition  period  is  about  19  hours. 

Suspended  below  the  glass  sphere  is  an  acoustic  transducer  with  source  strength  192  dB  re  1 
jiPa  at  1  meter  that  generates  and  receives  8-kHz  tone  bursts.  Pulses  10  ms  long  are  transmitted  by 
the  floats  in  a  programming  sequence.  A  different  float  transmits  every  45  seconds.  When  12 
floats  are  deployed,  each  float  transmits  every  9  minutes.  The  floats  are  listening  whenever  they 


14 


Figure  2.1  Schematic  drawing  of  a  typical  Swallow  float. 
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are  not  transmitting.  They  receive  pulses  transmitted  by  other  floats  as  well  as  surface/bottom  re¬ 
flections  of  their  own  pulses.  The  interfloat  and  fioat-to- surface  acoustic  travel  times  can  be  used 
to  localize  the  floats  as  well  as  to  correct  for  float  clock  drift. 


2.3  1989  Experiment  Summary 

The  12  Swallow  floats  were  deployed  for  a  24-hour  period  on  8  -  9  July  1989  near  34°  50'  N, 
122°  20'  W,  about  150  km  west,  northwest  of  Pt.  Arguello,  California.  The  Swallow  float  deploy¬ 
ment  was  a  part  of  the  Downslope  Conversion  experiment  that  has  as  its  primary  objective  the 
study  of  the  physics  of  downslope  signal  propagation. 

Of  the  12  floats,  9  were  freely  drifting  in  the  water  column,  and  3  were  tethered  to  the  ocean 
bottom  by  3.05-meter  lines  with  10-  to  15-  lb  anchors.  The  three  bottom-tethered  floats  were  posi¬ 
tioned  at  the  comers  of  a  triangle  with  sides  about  6.3  km  long  in  order  to  provide  an  absolute  ref¬ 
erence  for  the  float  localization.  The  nine  midwater  floats  were  deployed  in  a  quasi-vertical  line 
array  geometry  with  a  vertical  float  separation  of  about  400  m,  starting  at  about  600-m  depth  to 
about  3800  m.  The  midwater  floats  were  put  into  the  water  at  about  the  geometric  center  of  the 
bottom-float  triangle.  Figure  2.2  shows  the  Swallow  float  deployment  geometry  and  retrieval  po¬ 
sitions  of  all  12  Swallow  floats.  The  symbols  used  for  plotting  the  float  retrieval  positions  are  the 
float  identification  numbers.  Figures  2.3  shows  the  planned  float  deployment  depths  while  Table 
2.1  lists  the  planned  and  actual  deployment  depths  of  the  Swallow  floats.  In  this  deployment,  all 

Table  2.1  Planned  and  actual  deployment  depths,  July  1989  experiment. 


Boat 

Number 

Planned 
Deployment 
Depth  (m) 

Actual 
Achieved 
Depth  (m) 

0 

500 

560  -  580 

1 

900 

995-  1010 

2 

1300 

1375-  1395 

3 

1700 

1660-  1680 

4 

2100 

2080-2100 

5 

2500 

2495-2510 

6 

2900 

2915-2940 

7 

3300 

3350  -  3365 

8 

3700 

3770  -  3780 

9 

bottom 

4039 

bottom 

4051 

11 

bottom 

4041 

16 


Figure  2.2  Swal 
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Sound  Speed  (m/sec) 


Figure  2.3  Planned  Swallow  float  deployment  depth.  Also  plotted  in  this  figure  is  the 
sound  speed  profile  based  on  historical  temperature  and  salinity  data. 
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12  floats  that  were  taken  on  the  trip  were  deployed,  and  all  recorded  full  data  tapes  except  float  4, 
which  quit  recording  about  4  hours  earlier  than  expected  due  to  a  cassette  tape  defect.  Also,  float 
5’s  time  release  occurred  at  08:00, 9  July,  rather  than  at  08:00, 10  July,  due  to  a  programming  er¬ 
ror  that  caused  the  float  to  ascend  back  to  the  surface  prematurely.  In  all,  95  percent  of  the  poten¬ 
tial  number  of  records  that  could  have  been  recorded  by  12  properly  functioning  floats  were 
recorded. 

2.4  Data  Description 

2.4.1  8-KHz  Range  Data 

As  discussed  in  section  2.2,  the  floats  transmit  and  receive  10  msec,  8-tHz  pulses  in  order  to 
measure  interfloat  and  float-to-surface  acoustic  travel  times.  A  given  float  receives  pulses  trans¬ 
mitted  by  other  floats  as  well  as  its  own  arrivals  (surface  echoes,  bottom  bounces,  or  mixtures  of 
both).  Since  only  direct  path  pulses  between  floats  and  surface  echo  pulses  from  one  float  to  itself 
are  of  interest,  an  edge  detector  program  was  used  to  detect  and  extract  such  pulses  that  corre¬ 
spond  to  the  first  arrivals  in  a  narrow  range  time  window  [Hodgkiss  and  Anderson,  1983;  Culver 
andHodgkiss,  1988]. 

Figure  2.4  shows  the  leading  edge  of  float  l’s  detection  of  its  own  surface  echo  8-kHz  pulses 
and  pulses  transmitted  by  other  floats  according  to  its  own  clock  as  a  function  of  time  during  the 
July,  1989  experiment.  The  vertical  axis  is  converted  from  travel  times  in  seconds  to  range  (or 
depth)  in  meters,  i.e.,  the  direct  pulses  arrival  times  are  converted  to  range  by  using  a  speed  of 
sound  of  1500  m/s  and  the  self  surface  echo  returns  are  converted  to  depth  by  using  half  of  a 
speed  of  sound).  The  horizontal  axis  is  time  in  units  of  Swallow  float  records.  Each  Swallow  float 
record  equals  45  seconds.  The  horizontal  dotted  line  around  1000  m  indicates  the  approximate 
depths  at  which  float  1  stabilized. 

2.4.2  Acoustic  VLF  Data 

The  power  spectral  estimates  from  data  collected  by  float  l’s  hydrophone  during  records  1040 
to  1680  (00:28  and  08:28  PST,  9  July  1989)  are  presented  in  Figure  2.5  in  spectrogram  format. 
The  spectral  estimates  were  made  by  dividing  a  40.96-second  piece  of  data,  gotten  from  the  44 
seconds  of  data  in  each  record  after  skipping  the  first  3.04  seconds  in  order  to  reduce  tape  recorder 
contamination,  into  seven  512-point  segments  (each  segment  is  10.24  seconds  long)  with  a  50 
percent  overlap  between  segments.  The  segments  were  Fourier  transformed  after  being  windowed 
with  a  Kaiser- Bessel  window  of  a  -  2.5.  The  seven  spectra  were  then  incoherently  averaged  in  or¬ 
der  to  reduce  the  variance  of  the  power  spectral  estimates.  Prominent  features  of  this  data  set  in¬ 
clude 
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Figure  2.4  8-kHz  acoustic  range  pulses  received  by  float  1  (freely  drifting)  as  a  function 
of  time  during  the  July,  1989,  experiment.  The  horizontal  axis  of  the  figure  gives  the  time 
in  units  of  Swallow  float  record  number.  Eighty  records  represent  one  hour. 
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Record  Number 


Frequency  (Hz) 


Figure  2.5  VLF  acoustic  spectrogram  of  float  1, 00:28  -  08:28  PST,  9  July  1989.  The 
vertical  axis  of  the  figure  gives  the  time  in  units  of  Swallow  float  record  number. 

Eighty  records  represent  one  hour. 
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1.  MARK  VI  Source  Deployed  from  R/V  Aloha:  As  part  of  a  companion  experiment  being 
conducted  concurrently  with  the  Swallow  float  deployment,  a  MARK  VI  source  was 
towed  at  an  average  speed  of  8  knots  and  at  a  depth  of  about  90  meters  by  the  R/V  Aloha. 
This  source  ship  was  located  about  2500  km  and  10°  north  of  due  west  from  the  Swallow 
float  array  during  their  recording  period.  Figure  2.6  shows  the  MARK  VI  source  relative 
to  the  Swallow  float  deployment  site.  The  broadcast  plan  of  the  R/V  Aloha  was  to  transmit 
14  Hz  for  a  half-hour,  then  8  Hz  for  a  half-hour,  then  14  Hz,  then  1 1  Hz  and  then  repeat  the 
pattern.  Three  tonals  in  the  Swallow  float  frequency  band,  at  8  Hz,  11  Hz  and  14  Hz,  are 
visible  that  are  believed  to  be  generated  by  the  MARK  VI  source.  Evidently,  the  Swallow 
floats  can  see  quite  clearly  the  14-Hz  line  projected  by  the  MARK  VI  source  and,  at  times, 
the  1 1  Hz  and  the  8  Hz  lines. 

2.  Earthquake:  Around  record  1 535  (06:39, 9  July)  the  Swallow  float  detected  the  local  Rich¬ 
ter  magnitude  4.1  earthquake  that  occurred  along  the  Calaveras  fault  in  central  California. 

3.  Ship-Generated  Signals:  Harmonic  line  sets  broadcast  by  commercial  ships,  generated  by 
cavitation  resulting  from  the  low  pressure  of  the  props  that  provide  the  ship’s  thrust,  are  a 
well-known  contributor  to  the  infrasonic  sound  field  [Ross,  1976;  D’Spain  et.  al.,  1991  ]. 
The  Swallow  float  deployment  site,  just  west  of  the  mouth  of  the  Santa  Barbara  channel,  is 
an  area  of  very  heavy  shipping  traffic.  The  7-,  9.4-,  and  18-  to  19-Hz  spectral  lines  are  be¬ 
lieved  to  be  generated  by  surface  ships. 

2.4.3  Environmental  Data 

A  large  amount  of  environmental  data  was  collected  during  the  experiment  [Olivera,  1990], 
The  swell  was  observed  to  be  10  to  13  feet  in  height,  coming  from  330°  magnetic.  The  tempera¬ 
ture  was  in  the  low  sixties  in  Fahrenheit.  The  weather  in  this  part  of  the  ocean  comes  from  the 
northwest,  and  the  wind  speed  is  around  10  to  15  knots. 

Water  temperature  and  salinity  measurements,  collected  during  the  Downslope  Conversion 
and  a  companion  experiments,  relevant  to  this  study  include  (1 )  Twelve  expendable  bathythermo¬ 
graph  (XBT)  probes,  launched  near  the  Swallow  float  experiment  site  during  8  to  10  July,  reached 
a  depth  of  approximately  1800  meters;  (2)  One  conductivity -temperature-depth  measurement 
(CTD)  was  collected  on  1 1  July,  1989,  at  the  center  of  the  float-triangle  where  all  the  freely  drift¬ 
ing  floats  were  put  into  the  water.  The  CTD  station  sampled  from  surface  to  bottom  so  that  a  com¬ 
plete  sound  speed  profile  for  the  entire  water  column  could  be  obtained;  (3)  One  CTD 
measurement  to  the  depth  of  2000  meters  was  made  at  33°  48'  N,  141 0  03'  W  on  1 1  July,  1989; 
and  (4)  116  air-launched  bathythermographs  (AXBT),  designed  to  reach  a  depth  of  800  meters, 
were  dropped  by  a  P-3  aircraft  between  8  - 10  July,  1989,  along  the  path  between  32°  N,  150°  W 
and  35°  N,  122°  W  as  shown  in  Figure  2.6  marked  with  the  A’s.  Figure  2.7  is  the  complete  sum- 
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Figure  2.6  MARK  VI  source  -  Swallow  float  array  geometry,  9  July  1989.  The  line 
of  triangles  marks  the  AXBT  measurements  taken  between  8  - 10  July  1989. 
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mary  of  all  sound  speed  profiles  derived  from  the  measurements  taken  by  both  the  Downslope 
Conversion  and  a  companion  experiments  during  8-11  July  1989. 

2.5  Summary 

This  chapter  described  the  Swallow  float  system  and  summarized  the  Swallow  experiment 
conducted  on  8  -  9  July  1989.  MPL’s  12  Swallow  floats  were  deployed  for  a  24-hour  period  near 
34°  50'  N,  122°  20'  W,  approximately  150  km  west-northwest  of  Pt.  Arguello,  California.  Of  the 
12  floats,  9  were  freely  drifting  in  the  water  column  while  3  were  anchored  to  the  ocean  bottom  to 
provide  an  absolute  reference  for  the  float  navigation.  The  freely  drifting  floats  were  deployed  to 
span  the  water  column  of  4100-m  depth  with  a  vertical  float  separation  of  about  400  m,  starting  at 
about  600-m  depth  to  about  3800  m.  Data  collected  during  the  experiment  included  (1)  range 
data:  Floats  transmitted  and  received  8-kHz  pulses  in  a  preprogrammed  schedule  to  measure  in¬ 
terfloat  and  float-to- surface  travel  times  from  which  the  float  positions  could  be  determined  with  a 
postprocessing  float  localization  technique  based  on  the  least-squares  method;  (2)  VLF  acoustic 
data:  Floats  recorded  acoustic  data  in  the  1-  to  25-Hz  band.  Of  particular  interest  in  the  data  set 
was  a  14-Hz  tonal  generated  by  the  MARK  IV  source  deployed  by  the  R/V  Aloha,  a  source  ship 
involved  in  a  companion  experiment  located  about  2500  km  and  10°  north  of  due  west  from  the 
Swallow  float  array.  The  tonal  was  clearly  heard  by  the  floats  through  most  of  the  experiment;  and 
(3)  environment  data:  A  large  volume  of  environmental  data  such  as  AXBT,  XBT,  and  CTD  mea¬ 
surements  was  collected  between  34°  50'  N,  122°  20'  W  and  32°  00'  N,  150°  00'  W  during  8-10 
July  1989. 
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3.  SWALLOW  FLOAT  LOCALIZATION 


3.1  Introduction 

The  Swallow  float  array’s  geometry  must  be  determined  before  the  VLF  acoustic  data  can  be 
coherently  combined  to  form  a  beamformer.  As  described  in  last  chapter,  the  float  transmits  and 
receives  10-msec,  8-KHz  pulses  in  a  preprogrammed  sequence  to  measure  the  interfloat  and  sur¬ 
face  echo  travel  times.  Figure  3.1  is  a  conceptual  diagram  depicting  a  realization  of  the  geometry 


Sea  Surface 


Figure  3.1  Swallow  float  array  geometry  conceptual  diagram. 


of  the  freely  drifting  array.  The  arrow  lines  represent  the  travel  time  measurements  between  two 
floats  and  the  float-to-surface  echoes,  while  the  circles  represent  the  Swallow  floats.  The  focus  of 
this  chapter  is  to  estimate  the  Swallow  float  array  geometry  in  a  3-D  coordinate  system  from  the 
travel  time  measurements  as  a  function  of  time  during  the  July  1989  experiment. 

This  chapter  is  organized  as  follows.  Following  this  introductory  section,  Section  3.2  reviews 
the  theoretical  basis  for  float  localization  techniques.  The  desire  here  is  to  provide  the  necessary 
background  for  discussion  in  the  subsequent  sections.  Section  3.3  describes  the  inputs  to  the  float 
localization  filters,  while  Sections  3.4  and  3.5  discuss  approaches  to  filter  tuning  and  the  outputs 


26 


from  the  filters,  respective’y.  Lastly,  a  summary  of  the  float  localization  results  is  given  in  Section 
3.6.  Detailed  information  about  the  float  localization  is  documented  in  [Chen  and  Hodgkiss, 
1990]. 

3.2  Float  Localization  Techniques 

The  function  of  the  float  localization  filter  is  to  algorithmically  compute  float  positions  from 
noisy  range  measurements  in  a  statistical  sense  depending  on  the  optimality  criterion  chosen.  In 
this  section,  two  such  filters  developed  at  MPL  are  reviewed  from  the  performance  criteria  point 
of  view. 

3.2.1  Least-Squares  Filter 

The  generalized  (or  weighted)  least-squares  filter  [Gelb,  1974;  Sorenson,  1980]  involves  us¬ 
ing  the  current  set  of  range  measurements,  Zn  that  are  related  to  the  set  of  float  positions  XJ(  by 
the  expression 


Z  n  =  h(Xn)+Vn  (measurement  model)  (3.1) 

where  h  is  the  nonlinear  function  that  gives  the  ideal  (noiseless)  connection  between  the  range 
measurements  and  the  float  positions,  and  Vn  is  the  set  of  range  measurement  errors  that  are  also 
known  as  the  residuals.  The  least-squares  method  is  concerned  with  determining  the  most  proba¬ 
ble  set  of  float  positions  that  are  defined  as  the  set  of  float  positions  that  minimizes  the  weighted 
sum  of  the  squares  of  the  residuals: 


minimize  [Z„  -  *  (X,)  lX  [Z„  -  *  (X,)  ]  (3.2) 

X„ 

where  W;I  is  the  weighting  matrix  chosen  as  the  inverse  of  the  range  covariance  matrix  Rn,  de¬ 
fined  as 


E[VnV^]  =  (measurement  statistics)  (3.3) 

Thus,  the  least-squares  filter  estimates  current  float  positions  by  using  only  the  current  set  of 
measurements  and  an  estimate  of  the  measurement  error  statistics. 

3.2.2  Kalman  Filter 

The  Kalman  filter  [Gelb,  1974;  Brown,  1983;  Sorenson,  1985;  Culver  and  Hodgkiss,  1988]  at¬ 
tempts  to  estimate  float  positions  better  by  taking  advantage  of  the  knowledge  of  float  dynamics 
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and  incorporating  a  fading  memory  of  past  data  into  the  estimator  structure.  In  addition  to  the 
measurement  model  (3.1)  and  statistics  (3.3),  the  Kalman  filter  also  uses  the  system  model  (float 
dynamics),  the  system  statistics,  and  the  initial  conditions.  The  system  model  is  given  by 

Xn  =  <I>Xn_  j  +  rx„_  j  (system  model)  (3.4) 

where  <I>  is  the  state  transition  matrix  that  relates  Xn_  j  to  Xn,  r  is  the  matrix  that  relates  the  ac¬ 
celerations  to  the  float  positions,  and  X„  _  t  is  the  set  of  accelerations  with  zero  means  and  sec¬ 
ond-order  statistics  described  by 

E[X„Xrm]  =  Q„8nm  (system  statistics)  (3.5) 

The  initial  conditions  are 


So  “  E(X0) 

(3.6) 

P0=  e[(X0-*0)  (Xo-S,,)1) 

(3.7) 

where  X0  is  the  estimate  of  the  initial  float  positions,  and  P0  is  the  estimate  of  the  error  variance 
in  the  initial  float  positions. 

The  operation  of  the  Kalman  filter  can  be  viewed  as  a  predictor-corrector  process  [Candy, 
1988].  First,  suppose  we  have  available  to  us  the  previous  float  position  estimates  Xn  _ ,  and  asso¬ 
ciated  position  covariance  matrix  Pn  _  x : 

P„-l  =  e[(X„_,-X„_1)  (X„_. (3.8) 

and  we  would  like  to  obtain  the  best  estimate  of  the  current  float  positions  based  on  the  previous 
estimate.  We  are  in  the  “prediction  phase”  of  the  process.  The  system  model  (3.4)  is  used  to  pre¬ 
dict  the  float  positions  and  the  associated  position  covariance  matrix: 

*„(-)«**„_,  (3.9) 

p„(-)  =  &p„_)<pT+rQrr  (3.io) 

The  “minus”  is  introduced  here  (following  the  notation  in  reference  [Gelb,  1974])  as  a  reminder 
that  this  is  our  best  estimate  prior  to  incorporating  the  measurement.  We  now  seek  to  use  the  mea¬ 
surement  Zn  to  improve  the  estimate  X„(-).  We  choose  a  linear  blending  of  the  noisy  measure¬ 
ment  and  the  %.„(-)  in  accordance  with  the  equation: 
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(3.11) 


X„  =  X„(-)  +  Kn  [Zn  -  h  (X„(-) )  ] 

The  optimal  choice  of  Kn,  the  Kalman  gain,  is  to  minimize  a  scalar  sum  of  the  diagonal  elements 
of  the  position  covariance  matrix  P„: 


minimize  {trace  [ PJ} 

K 


(3.12) 


where 

P„  =  E[(X,-X,.l)(X,-f(,.l)’]  (3.13) 

K;I  is  obtained  as 

K„  -  P,«H,r[H,P,(-)H„r+R„]"'  (3.14) 

or  in  a  simpler  form 

K„  =  P.HX'  (3-15) 

where 


H 


n 


(3.16) 


We  then  use  the  current  set  of  range  measurements  Zn  to  determine  the  innovations.  The  innova¬ 
tion  is  defined  as  Zn  -  h  [X„(-)  ] ,  which  is  the  difference  between  the  actual  range  measurement 
and  the  predicted  range  measurement  (3.9).  Now  we  enter  the  “correction  phase”  of  the  process. 
That  is,  we  correct  or  update  the  predicted  positions  based  on  the  new  information  in  the  range 
measurements  -  the  innovation.  The  innovation  is  weighted  by  the  Kalman  gain  Kn  to  correct  the 
predicted  positions  &„(-)  as  described  in  (3.11).  The  associated  position  covariance  matrix  (3.10) 
is  corrected  as  well: 


P„  =  n-K„H„]P„(-)  (3.17) 

The  predictor-corrector  process  then  repeats  until  all  measurements  are  consumed. 

A  useful  interpretation  [Gelb,  1974]  of  the  Kalman  gain  (3.15)  is  that  the  gain  is  “proportion¬ 
al”  to  the  uncertainty  in  the  position  estimates,  and  “inversely  proportional”  to  the  range  measure¬ 
ment  noise.  In  other  words,  if  the  range  measurement  errors  are  large  and  the  predicted  position 
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errors  are  small,  the  innovation  is  due  chiefly  to  the  noise,  and  only  a  small  change  in  the  predict¬ 
ed  positions  should  be  made.  On  the  other  hand,  small  measurement  noise  and  large  uncertainty  in 
the  position  estimates  suggest  that  the  innovation  contains  considerable  information  about  errors 
in  the  position  estimates.  Therefore,  the  difference  between  the  actual  and  the  predicted  range 
measurement  will  be  used  as  a  basis  for  strong  corrections  to  the  position  estimates. 

In  short,  the  Kalman  filter  estimates  float  positions  by  weighting  a  current  set  of  measure¬ 
ments  against  previous  position  estimates  propagated  forward  in  time  using  an  equation  of  mo¬ 
tion.  Thus,  the  Kalman  filter  uses  all  of  the  following:  current  measurements  as  well  as  past 
measurements  in  a  fading  memory  fashion,  estimates  of  the  measurement  error  statistics  and  ac¬ 
celeration  error  (process  noise)  statistics,  and  the  float  dynamics.  It  has  been  shown  [Culver  and 
Hodgkiss,  1988]  that  the  Kalman  filter  outperforms  the  least-squares  filter  in  the  presence  of  high 
measurement  noise  and  when  the  process  noise  is  low  because  of  its  ability  to  track  float  motion 
and  effectively  smooth  noisy  measurements. 

3.3  Inputs  to  the  Localization  Filter 

This  section  presents  the  inputs  to  both  filters  using  data  from  the  July  1989  experiment. 
Based  on  the  discussion  in  the  previous  section,  input  data  fall  under  four  categories:  measure¬ 
ment  data,  measurement  statistics,  initial  estimates,  and  process  noise. 

3.3.1  Measurement  Data 

Range  measurements  used  by  the  filters  are  derived  quantities  computed  from  the  8-kHz  pulse 
travel  time  measurements  between  floats  (or  from  float  to  surface)  using  an  estimate  of  the  sound 
speed.  Thus,  measurement  data  include  travel  time  estimates  and  sound  speed  estimates. 

3.3.1. 1  Travel  Time  Estimates 

Measurement  of  pulse  travel  times  has  received  a  great  deal  of  attention  during  the  process  of 
the  Swallow  float  system  development  in  the  past  few  years.  References  [Hodgkiss  and  Anderson, 
1983;  Culver  et.  al.,  1989]  have  addressed  this  subject  in  detail.  In  brief,  each  float  is  equipped 
with  an  acoustic  transducer  that  transmits  and  receives  8-kHz  pulses.  Pulses  are  transmitted  by  the 
floats  in  a  preprogrammed  sequence.  A  different  float  transmits  every  45  seconds;  12  floats  were 
deployed  in  this  experiment,  and  each  float  transmits  every  9  minutes.  A  given  float  receives  puls¬ 
es  transmitted  by  other  floats  as  well  as  its  own  arrivals  (surface  echoes,  bottom  bounces,  or  mix¬ 
tures  of  both).  Since  only  direct  path  pulses  between  floats  and  surface  echo  pulses  from  one  float 
to  itself  are  of  interest,  an  edge  detector  program  was  used  to  detect  and  extract  such  pulses  that 
correspond  to  the  first  arrivals  in  a  narrow  range  time  window. 
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Surface  echo  travel  time  measurements  are  obtained  by  subtracting  outgoing  pulse  transmit 
times  from  pulse  arrival  times,  whereas  interfloat  travel  time  measurements  are  obtained  by  sub¬ 
tracting  pulse  transmit  time  according  to  the  transmitting  float  from  pulse  arrival  times  at  the  re¬ 
ceiving  float.  Since  two  time  bases  are  involved,  interfloat  travel  time  measurements  contain  a 
bias  due  to  variation  in  the  float  clock  rates.  It  has  been  shown  [Culver  and  Hodgkiss,  1989]  that 
the  bias  can  be  reduced  significantly  by  averaging  reciprocal  path  interfloat  travel  time  measure¬ 
ments.  Because  interfloat  measurements  were  not  made  simultaneously,  travel  time  measurements 
must  be  interpolated  by  a  factor  of  12  before  reciprocal  path  measurements  may  be  averaged. 

Figures  3.2  and  3.3  show  the  typical  direct  path  pulse  detections  between  two  freely  drifting 
floats  (2  and  8)  whereas  Figure  3.4  shows  the  interfloat  travel  time  estimates  calculated  by  averag¬ 
ing  interpolated  reciprocal  path  measurements.  Figure  3.5  shows  the  leading  edge  of  each  float’s 
detection  of  its  own  transmitted  and  surface  echo  8-kHz  pulses.  The  vertical  axis  in  the  figure  has 
been  scaled  from  travel  time  to  depth,  using  half  of  1500  m/s  for  sound  speed.  The  lines  of  dots 
correspond  to  each  float’s  detection  of  the  incoming  pulses,  i.e.,  surface  echo  pulses.  Note  in  Fig¬ 
ure  3.5  that  the  lines  of  dots  for  floats  4  and  5  terminate  prematurely.  This  is  due  to  the  fact  that 
float  4  stopped  recording  data  4  hours  too  early,  and  float  5’s  automatic  ballast  release  time  was 
set  incorrectly  (as  indicated  in  section  2.3).  Figure  3.6  is  the  enlarged  version  of  the  surface  echo 
travel  time  measurements  for  the  bottom  float  9.  As  pointed  out  in  [Culver  and  Hodgkiss,  1989], 
the  bottomed  floats’  surface  measurements  like  these  are  expected  to  be  approximately  constant 
because  the  floats  themselves  are  approximately  stationary;  they  are  tethered  to  the  ocean  bottom 
by  3.05-meter  lines.  Noise  in  the  surface  echo  pulse  arrival  time  is  thought  to  be  caused  by  scat¬ 
tering  of  the  pulse  at  the  rough,  moving  sea  surface  and  destructive  interference  among  multiple 
arrivals  at  the  receiver.  Since  filter  performance  can  be  improved  by  substituting  a  constant  for  the 
actual  measurements,  the  most  probable  value,  i.e.,  the  mode  of  the  travel  time  measurements  as 
shown  in  Figure  3.7  (5.402  seconds  for  float  9;  5.422  seconds  for  float  10;  and  5.404  seconds  for 
float  11)  will  be  used  as  the  travel  time  estimates.  Table  3.1  lists  the  modes  of  the  surface  echo 
travel  time  estimates  for  the  bottomed  floats. 

Table  3.1  Surface  echo  travel  time  estimates  for  bottomed  floats. 


Float 

Number 

Surface  Echo 
Estimate 
(msec) 

9 

5402 

10 

5422 

11 

5404 
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Figure  3.2  Float  2  listening  to  float  8,  July  1989  experiment.  The  horizontal  axis  of 
the  figure  gives  the  time  in  units  of  Swallow  float  record  number.  Eighty  records 

represent  one  hour. 
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Figure  3.3  Float  8  listening  to  float  2,  July  1989  experiment.  The  horizontal  axis  of 
the  figure  gives  the  time  in  units  of  Swallow  float  record  number.  Eighty  records 

represent  one  hour. 
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Figure  3.4  Travel  time  estimates  between  two  freely  drifting  floats  (2  and  8).  The 
horizontal  axis  of  the  figure  gives  the  time  in  units  of  Swallow  float  record  number. 

Eighty  records  represent  one  hour. 
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Figure  3.5  Leading  edges  of  each  float’s  detection  of  its  own  surface  echo  pulses.  The 
horizontal  axis  of  the  figure  gives  the  time  in  units  of  Swallow  float  record  number. 

Eighty  records  represent  one  hour. 
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Figure  3.6  Float  1 1  surface  echo  pulses,  July  1989  experiment.  The  horizontal  axis  of 
the  figure  gives  the  time  in  units  of  Swallow  float  record  number.  Eighty  records 

represent  one  hour. 


Figure  3.7  Histogram  of  float  1 1  surface  echoes  minus  the  mean  541 1  msec  (binwidth 

is  one  msec),  July  1989  experiment. 


Figures  3.8  and  3.9  show  the  surface  reflected  pulse  detections  between  the  bottomed  floats,  9 
and  1 1 .  The  direct  path  pulses  between  the  bottomed  floats  were  not  detected  because  the  sound 
speed  gradient  bends  nearly  all  of  the  sound  energy  upward  over  the  6.3-km  distance  that  sepa¬ 
rates  the  floats.  Direct  path  travel  times  can  be  calculated  from  surface  reflection  travel  times 
[Culver  and  Hodgkiss,  1989]: 
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direct  path 
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(3.18) 


where  t”  is  the  interfloat  travel  time,  x.(  and  t..  are  the  individual  float  surface  echo  travel  times. 
chmss *s  harmonic  mean  sound  speed  [Polvani,  1984]  for  a  vertical  path  from  the  surface  to  the 
bottom,  and  c,--  is  the  sound  speed  at  the  depth  of  floats  i  and  j. 


Like  the  bottom  float  surface  echo  travel  times,  the  interfloat  travel  time  between  two  bottom 
floats  should  be  approximately  constant  because  the  floats  are  approximately  stationary.  Figure 
3.10  shows  the  direct  path  travel  time  estimates,  calculated  using  equation  (3.18)  and  the  constant 
value  (mode)  estimates  of  the  averaged  interpolated  reciprocal  surface  reflection  travel  times  be¬ 
tween  the  bottom  floats  9  and  10.  Because  the  underlying  direct  path  travel  time  must  be  approxi¬ 
mately  constant,  the  modes  of  the  direct  path  travel  time  estimates  will  be  used  the  as  the  direct 
path  travel  times  estimates  as  shown  in  Figure  3.11.  Table  3.2  lists  the  modes  of  the  direct  path 
travel  times  between  the  bottomed  float  pairs. 


Table  3.2  Direct  path  travel  time  estimates  between  bottomed  floats. 


Float  Path 

Direct  Path 
Estimate 
(msec) 

9-10 

4025 

9-11 

4157 

10-11 

4132 

3.3.1.2  Sound  Speed  Estimates 

The  sound  speed  profile  for  the  upper  1000  meters  was  derived  from  five  expendable 
bathythermograph  (XBT)  measurements  made  from  the  R/V  New  Horizon  near  the  deployment 
site.  The  measurements  were  made  at  various  times  during  the  Swallow  float  experiment  and  his¬ 
torical  salinity  data  archived  by  the  National  Oceanographic  Data  Center  [Churgin  and  Halmins- 
ki,  1974]  for  area  24,  3rd  quarter  of  the  year,  in  the  upper  1000  meters  were  used  in  deriving  the 
sound  speed  profiles.  The  sound  speed  profile  for  the  lower  3 1 00  meters  was  obtained  from  a  con- 
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Figure  3.8  Float  9  listening  to  float  1 1,  July  1989  experiment.  The  horizontal  axis  of 
the  figure  gives  the  time  in  units  of  Swallow  float  record  number.  Eighty  records 

represent  one  hour. 
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Figure  3.9  Float  11  listening  to  float  9.  July  1989  experiment.  The  horizontal  axis  of 
the  figure  gives  the  time  in  units  of  Swallow  float  record  number.  Eighty  records 

represent  one  hour. 
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Figure  3.10  Direct  path  travel  time  estimates  between  bottomed  floats  9  and  1 1 , 
calculated  using  constant  mode  values  of  surface  echo  travel  times,  July  1989 
experiment.  The  horizontal  axis  of  the  figure  gives  the  time  in  units  of  Swallow 
float  record  number.  Eighty  records  represent  one  hour. 


Figure  3.11  Histogram  of  direct  path  travel  time  estimates  minus  the  mean  41 59  msec 
(bin width  is  one  msec),  July  1989  experiment. 
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ductivity-temperature-depth  CTD  station  [Olivera,  1990]  carried  out  on  July  11,  1989,  near  the 
center  of  the  float-triangle  where  all  the  freely  drifting  floats  were  put  into  the  water.  Using  the 
equation  given  by  Mackenzie  [Mackenzie,  1981],  the  composite  sound  speed  profile  was  calculat¬ 
ed.  Harmonic  mean  sound  speeds  for  each  interfloat  and  float  to  surface  path  were  also  calculated 
and  are  given  in  Appendix  A. 

3.3.2  Measurement  Statistics 

The  measurement  error  variance  estimates  are  used  to  weight  the  measurements  so  that  great¬ 
er  weight  is  given  to  measurements  with  a  smaller  error  variance.  Since  travel  time  measurement 
error  and  sound  speed  estimate  error  are  assumed  to  be  uncorrelated,  the  measurement  error  vari¬ 
ance  for  float  pair  ij  can  be  expressed  as  [Culver,  1988]: 


where  t".  is  the  estimated  travel  time  between  float  i  and  j,  o2  is  the  sound  speed  error  variance 
in  c. Cjj  is  the  estimated  harmonic  sound  speed  between  float  i  and  j,  and  o2  is  the  travel  time 

J  J  Cjj 

error  variance  between  float  i  and  j. 

3.3.2.1  Travel  Time  Variances 

It  has  been  shown  [Culver  et.  al.,  1989]  that  the  variance  of  a  bottomed  float’s  surface  echo 
travel  time  measurements  can  be  used  as  an  estimate  of  the  variance  of  the  measurement  error,  be¬ 
cause  the  true  surface  echo  times  are  approximately  constant.  Subtracting  the  means  from  the  sur¬ 
face  echo  travel  time  estimates  produces  estimates  of  the  travel  time  estimate  errors.  The 
variances  of  the  error  time  series  are  given  in  Table  3.3. 


Table  3.3  Estimated  variance  of  bottomed  float  travel  time,  July  1989  experiment. 


Float 

Number 

Float 

Number 

Variance 

(msec2) 

9 

9 

115 

10 

10 

110 

11 

11 

97 

9 

10 

78 

9 

11 

64 

10 

11 

142 
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In  the  previous  section,  the  bottomed  float  surface  echo  travel  time  measurements  were  ruled 
too  noisy  to  be  used  directly  in  estimating  the  float  positions,  and  the  mode  of  the  travel  time  mea¬ 
surements  was  taken  as  the  travel  time  estimate.  The  variance  of  the  error  in  the  constant  estimate 
(mode)  is  expected  to  be  much  smaller  than  the  variance  of  the  error  in  the  measurement.  Based 
upon  the  predicted  maximum  vertical  movement  of  a  tethered  bottom  float  and  wave  height  in  the 
sea  surface  during  the  deployment  (about  3  meters),  the  error  in  the  mode  is  taken  to  be  mean  zero 
and  to  have  a  variance  of  6  msec2.  By  the  same  token,  the  interfloat  travel  time  variance  is  taken 
to  be  4  msec2  [Culver  and  Hodgkiss,  1989]  based  upon  the  predicted  maximum  horizontal  move¬ 
ment  of  a  bottom  float.  Table  3.4  summarizes  the  estimated  error  variances  of  the  bottomed  float 
travel  time  estimates. 


Table  3.4  Estimated  variance  of  bottomed  float  travel  time  based  on  predicted  float  movement, 

July  1989  experiment. 


Float 

Number 

Float 

Number 

Variance. 

msec2 

9 

9 

6 

10 

10 

6 

11 

11 

6 

9 

10 

4 

9 

11 

4 

10 

11 

4 

The  surface  echo  travel  time  estimate  error  for  the  freely  drifting  floats  cannot  be  estimated  di¬ 
rectly.  Because  these  floats  were  deployed  shallower  than  the  bottom  floats,  their  variances  are  ex¬ 
pected  to  be  lower  than  the  values  in  Table  3.3  but  higher  than  those  in  Table  3.4.  An  arbitrary,  yet 
logical  choice  of  error  variances  for  the  freely  drifting  float  surface  echo  travel  time  estimates  is 
listed  in  Table  3.5. 

Travel  time  error  variance  for  float  pairs  in  which  one  of  the  floats  is  not  stationary  can  be  es¬ 
timated  from  the  difference  between  reciprocal  path  travel  time  measurements  [Culver  et.  al., 
1989].  The  difference  in  clock  rates  causes  the  travel  time  difference  to  be  a  linear  function  of 
time.  The  difference  also  appears  to  contain  second-order  time  dependence  due  to  the  difference 
in  clock  accelerations.  Subtracting  a  second-order  fit  from  the  travel  time  difference  produces  an 
estimate  of  the  error  in  the  travel  time  estimates.  Figures  3.12  and  3.13  contain  the  interfloat  travel 
time  differences  for  floats  2  and  8  along  with  second-order  curves  fitted  to  the  differences  and  the 
corresponding  error  estimate  time  series.  The  variances  of  the  error  time  series  for  all  float  pairs 
are  summarized  in  Appendix  B. 
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Figure  3.12  Travel  time  difference  between  floats  2  and  8  and  a  second-order  fit,  July 
1989  experiment.  The  horizontal  axis  of  the  figure  gives  the  time  in  units  of  Swallow 
float  record  number.  Eighty  records  represent  one  hour. 


Table  3.5  Estimated  variance  of  freely  drifting  float  travel  time,  July  1989  experiment. 


Float 

Numbers 

Float 

Number 

Variance, 

msec2 

0 

0 

6 

1 

1 

10 

2 

2 

20 

3 

3 

30 

4 

4 

40 

5 

5 

50 

6 

6 

60 

7 

7 

70 

8 

8 

80 

3.3.2.2  Sound  Speed  Variances 

It  has  been  shown  [Culver,  1988;  Culver  and  Hodgkiss,  1989]  that  the  variance  of  the  sound 
speed  estimate  at  a  particular  depth  is  just  under  0.4  (meter/ sec)2  in  the  North  Pacific  Ocean.  The 
variance  of  the  harmonic  mean  sound  speed  for  a  given  path  would  be  expected  to  be  less  than 
that  in  the  sound  speed  at  a  particular  depth  due  to  the  averaging  effect.  The  error  in  the  harmonic 
mean  sound  speed  is  estimated  to  be  0.1  (meter/sec)2  for  the  1989  experiment. 

3.3.3  Initial  Estimates 

The  MPL  implementation  of  the  filters  requires  an  initial  estimate  of  the  float  positions.  Addi¬ 
tionally,  the  Kalman  filter  requires  an  initial  estimate  of  the  float  velocities  and  the  initial  position 
and  velocity  variance  estimates.  Before  estimating  the  initial  positions,  we  must  first  define  a  co¬ 
ordinate  system.  The  coordinate  system  for  the  float  localization  is  established  in  which  the  origin 
lies  at  the  surface  directly  above  float  9.  The  Z  axis  is  vertical,  positive  downwards  and  extends 
through  float  9,  and  the  X  axis  is  oriented  such  that  float  10  lies  in  the  XZ  plane.  The  positions  of 
floats  9, 10,  and  1 1  are  taken  to  be  the  ship’s  position  when  the  floats  were  launched  into  the  wa¬ 
ter.  Figure  3.14  shows  a  plane  view  of  the  coordinate  system.  The  locations  of  the  freely  drifting 
floats  will  be  estimated  relative  to  these  axes. 
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Figure  3.14  Plane  view  of  coordinate  system  used  for  float  localization. 


3.3.3.1  Initial  Position  Estimates 

The  time  for  the  initial  float  position  chosen  for  float  localization  is  record  1003,  00:00  July  9 
1989,  at  which  seven  of  the  nine  freely  drifting  floats  had  reached  equilibrium  depth.  The  MPL 
implementation  of  the  least-squares  filter  requires  a  rough  initial-position  estimate  to  bootstrap 
the  filter.  Using  depths  taken  from  Figure  3.5,  we  estimated  the  rough  initial  positions  to  be  the 
ship’s  position  relative  to  the  coordinate  system  when  the  floats  were  launched  into  the  water.  Ta¬ 
ble  3.6  lists  the  rough  initial-position  estimate. 
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Table  3.6  Rough  initial-position  estimate,  record  1003, 00:00  July  9  1989. 


Float 

Number 

X 

(meters) 

Y 

(meters) 

Z 

(meters) 

9 

0 

0 

4050 

6150 

0 

4050 

11 

3075 

5500 

4050 

3075 

2750 

570 

1 

3075 

2750 

990 

2 

3075 

2750 

1380 

3 

3075 

2750 

1670 

4 

3075 

2750 

2080 

5 

3075 

2750 

2490 

6 

3075 

2750 

2840 

7 

3075 

2750 

3350 

8 

3075 

2750 

3650 

Using  the  rough-initial  position  estimates  in  Table  3.6;  the  travel  time  estimates  for  record 
1003;  the  travel  time  variance  estimates  given  in  Tables  3.4,  and  3.5,  and  Appendix  B;  and  the 
sound  speeds  and  sound  speed  variances  given  in  Appendix  A,  the  least-squares  filter  produced 
the  position  estimate  given  in  Table  3.7.  It  has  a  root  mean  squared  (rms)  residual  of  2.67  meters. 
Since  an  initial-position  estimate  is  determined  to  be  satisfactory  when  the  least-squares  filter 
converges  to  a  position  estimate  that  results  in  an  rms  residual  of  less  than  7.5  meters  [Culver  and 
Hodgkiss,  1989],  the  position  estimate  at  record  1003  is  taken  as  a  good  initial-position  estimate 
and  will  be  used  by  the  Kalman  filter. 

3.3.3.2  Initial  Velocity  Estimates 

Using  the  good  initial-position  estimates,  the  least-squares  filter  was  again  run  for  the  first  13 
records  (9-minute  period)  from  records  1003  to  1015.  The  position  estimates  produced  by  the  fil¬ 
ter  were  used  to  calculate  the  position  change  rates,  i.e.,  (X,,  + ,  -  Xn)  /  45  metcrs/second. 
Twelve  such  successive  position  change  rates  were  then  averaged  to  produce  the  initial  velocity 
estimates.  Table  3.8  lists  the  velocity  estimates  for  record  1003, 00:00, 9  July  1989. 
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Table  3.7  Initial-position  estimate  produced  by  the  least-squares  filter  for  record  1003, 00:00, 

July  9  1989. 


Float 

Number 


9 

10 
11 


X 

(meters) 

Y 

(meters) 

Z 

(meters) 

0 

0 

4039 

6131 

0 

4051 

3103 

5518 

4041 

3779 

-1456 

574 

4887 

82 

999 

3885 

162 

1380 

4425 

969 

1674 

4110 

1398 

2081 

4092 

1482 

2489 

3870 

1754 

2841 

2875 

2085 

3353 

2824 

1871 

3650 

Table  3.8  Initial  velocity  estimate  (in  meters/second),  record  1003, 00:00,  July  1989  experiment. 
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3.3.3.3  Initial  Position  and  Velocity  Variance  Estimates 

The  estimate  of  the  variances  in  the  initial  positions  and  velocities  is  also  needed  by  the  Kal¬ 
man  filter.  The  estimate  suggested  by  [Culver  and  Hodgkiss,  1989]  will  be  used  and  is  listed  in 
Table  3.9. 

The  last  piece  of  information  required  by  the  Kalman  filter  is  the  float  acceleration  variance, 
also  known  as  the  process  noise,  described  in  equation  (3.5).  The  process  noise  is  a  parameter 
used  by  the  Kalman  filter  to  determine  how  much  weight  to  give  to  its  own  track  of  the  float  rela¬ 
tive  to  the  measurements.  When  the  model  (i.e.,  filter)  process  noise  matches  the  true  process 
noise  experienced  by  the  float  during  the  experiment,  the  Kalman  filter  is  performing  optimally 
and  will  have  more  estimates  falling  inside  the  predefined  confidence  interval  [Bar-Shalom  and 
Fortman,  1988],  However,  there  is  no  easy  way  of  estimating  the  true  process  noise.  One  ap¬ 
proach  suggested  by  [Culver  and  Hodgkiss,  1989]  is  to  run  the  filter  several  times,  each  time  with 
a  different  process  noise  value,  and  the  one  that  minimizes  the  innovation  power  will  be  selected 
as  the  candidate.  We  will  use  this  matched  processing  technique  with  values  ranging  from  1(T12 
to  5  x  1(T8  (meter/second2)2  in  search  of  the  true  process  noise. 


Table  3.9  Estimate  of  initial  position  and  velocity  variance  used  by  the  Kalman  filter,  July  1989 

experiment. 


Float 

Number 

Pos. 

Var. 

X 

(m2) 

Pos. 

Var. 

Y 

(m2) 

Pos. 

Var. 

z 

(m2) 

Vel.  Var. 
X 

(m/sec)2 

Vel. 

Var. 

Y 

(m/ 

sec)2 

Vel. 

Var. 

z 

(m/ 

sec)2 

9 

0.0 

0.0 

0.01 

0.0 

0.0 

0.0 

10 

0.01 

0.0 

0.01 

0.0 

0.0 

0.0 

11 

0.01 

0.01 

0.01 

0.0 

0.0 

0.0 

0 

0.01 

0.01 

0.01 

0.0001 

0.0001 

0.00001 

1 

0.01 

0.01 

0.01 

0.0001 

0.0001 

0.00001 

2 

0.01 

0.01 

0.01 

0.0001 

0.0001 

0.00001 

3 

0.01 

0.01 

0.01 

0.0001 

0.0001 

0.00001 

4 

0.01 

0.01 

0.01 

0.0001 

0.0001 

0.00001 

5 

0.01 

0.01 

0.01 

0.0001 

0.0001 

0.00001 

6 

0.01 

0.01 

0.01 

0.0001 

0.0001 

0.00001 

7 

0.01 

0.01 

0.01 

0.0001 

0.0001 

0.00001 

8 

0.01 

0.01 

0.01 

0.0001 

0.0001 

0.00001 
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3.4  Filter  Tuning 

3.4.1  Least-Squares  Filter 

Using  the  good  initial  positions  given  in  Table  3.7,  the  sound  speeds  and  sound  speed  varianc¬ 
es  given  in  Appendix  1,  the  travel  time  estimates  for  records  1003  to  2120,  and  the  travel  time 
variances  given  in  Tables  3.4  and  3.5  and  Appendix  2,  the  least-squares  filter  produced  position 
estimates  with  an  average  rms  of  2.7  meters.  The  residual  is  low  and  indicates  that  the  position  es¬ 
timate  is  acceptable.  However,  a  close  inspection  of  the  residual  sequences  [Zn  -  h  (\n)  ]  re¬ 
veals  that 


E  [Zd: 


direct  path  between  9  &  1 0 


-h  (4 


direct  path  between  9  &  10 


)]  = 


E[ 


afloat  9  surface  echo 


-A  (4 


oat  9  surface  echo 


)]  = 


1.3  meters 

(3.20) 

3.0  meters 

(3.21) 

1.2  meters 

(3.22) 

1.5  meters 

(3.23) 

4.9  meters 

(3.24) 

0.6  meters 

(3.25) 

By  definition,  residuals  in  the  measurement  model  are  to  be  mean  zero: 

E[VJ  =  E  [Zn  -  h  (Xn)  ]  =  0 


(3.26) 


The  bias  contained  in  the  residual  sequences  (3.20)  -  (3.25)  is  an  indication  that  the  mode  values 
selected  as  the  travel  time  estimates  may  not  be  the  optimal  choice  for  the  July  1989  data  set.  In 
order  to  improve  the  filter  performance,  trends  in  the  residual  sequences  must  be  removed.  An  it¬ 
erative  procedure  was  developed  to  achieve  this  and  is  described  as  follows: 

a-  _  E  by  -  h  (x”j)  ] 

1.  Compute  where  i -  9, 10, 11  and  j  «  9,  10, 11. 

2.  Update  =  t;.  -  At_ 

3.  Rerun  the  least-squares  filter. 

4.  If  all  |e  [z^-  h  {x ".)  ]  |  <  1  meter  then  stop;  otherwise  go  to  step  1 . 
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The  procedure  was  applied  to  the  July  1989  data  set.  After  10  iterations,  all  6  expected  residuals 
converged  to  within  1  meter,  and  the  least-squares  filter  produced  the  position  estimate  with  an 
average  rms  of  2.3  meters.  The  residual  is  lower  than  2.7  meters  and  indicates  that  the  position  es¬ 
timate  for  records  1003  to  2120  is  closer  to  the  true  float  positions.  The  constant  value  estimates 
used  in  the  10th  run  listed  in  Table  3.10  will  replace  the  mode  values  as  the  new  travel  time  esti¬ 
mates  and  will  be  used  by  the  Kalman  filter. 


Table  3.10  Bottomed  float  surface  echo  and  direct  path  between  bottomed  floats  travel  time 

estimate,  July  1989  experiment. 


Float-Path 

Travel  Time  Estimate 
(msec) 

9-10 

4028  (was  4025) 

9-11 

4152  (was  4157) 

10-11 

4130  (was  4132) 

9-9 

5401  (was  5402) 

10-10 

5409  (was  5422) 

11-11 

5405  (was  5404) 

3.4.2  Kalman  Filter 

The  measurements,  sound  speed  and  variance  estimates,  and  the  initial  estimates  were  applied 
to  the  Kalman  filter.  Using  process  noise  values  ranging  from  10~12  to  5  x  10~8  (meters/second)2, 
the  Kalman  filter  was  run  24  times.  Table  3.11  lists  the  mean  innovation  power  resulting  from 
each  run. 

It  has  been  shown  using  simulations  [Culver  and  Hodgkiss,  1989]  that  the  process  noise  that 
minimizes  the  power  in  the  innovations  sequence  produces  a  position  estimate  with  the  smallest 
true  rms  error.  The  minimum  is  seen  to  occur  at  1(T9,  which  is  the  same  found  for  the  1987  exper¬ 
iment  [Culver  and  Hodgkiss,  1989].  Thus,  process  noise  for  the  July  1 989  experiment  is  estimated 
to  be  1  O'9  (m/s2)2. 
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Table  3.11  Mean  innovation  power  produced  by  the  Kalman  filter,  July  1989  experiment 


Process  Noise 
((meters/ 
second2)2) 

Mean 

Innovation 

Magnitude 

(meters) 

1  X  I0~12 

7.0671 

2  X  10-12 

6.0969 

3  X  10~12 

5.6302 

5  X  10-12 

5.1290 

7  X  10“12 

4.8450 

1  X  10~U 

4.5793 

2  X  10-11 

4.1538 

3  X  10~U 

3.9536 

5  X  10“ 11 

3.7460 

7  X  10-11 

3.6335 

1  X  10~10 

3.5331 

2  X  10~10 

3.3866 

3  X  10" 10 

3.3269 

5  X  10"10 

3.2753 

7  X  10“ 10 

3.2546 

1  X  10-9 

3.2435 

2  X  10“9 

3.2539 

3  X  10“9 

3.2808 

5  X  10-9 

3.3401 

7  X  10-9 

3.3974 

1  X  to-8 

3.4770 

2  X  10-8 

3.6933 

3  X  10-8 

3.9648 

5  X  10-8 

6.8216 

3.5  Localization  Results 

In  this  section,  the  outputs  from  the  localization  filters  are  discussed.  Of  interest  are  RMS  po¬ 
sition  error  estimate,  RMS  residual/innovation,  float  position  estimate,  distance  difference  be¬ 
tween  filters,  float  depth  variation,  and  float  speed  estimate. 

The  curves  shown  in  Figure  3.15  are  the  mean  innovation  power,  or  the  Kalman  filter’s  esti¬ 
mate  of  the  rms  position  error,  and  the  least-squares  filter’s  estimate  of  the  rms  position  error.  The 
vertical  axis  is  rms  error  or  mean  innovation  magnitude  in  meters.  The  horizontal  axis  is  the  esti¬ 
mate  of  process  noise  given  to  the  Kalman  filter  as  a  parameter.  Simulations  [Culver,  1988]  have 
shown  that  the  time  rms  error  in  the  Kalman  filter’s  position  estimate  was  always  bracketed  by  the 
mean  magnitude  of  the  innovation  and  the  filter’s  estimate  of  the  rms  error.  Therefore,  the  rms  er¬ 
ror  in  the  Kalman  filter  estimate  is  thought  to  be  less  than  3.3  meters.  The  least-squares  filter’s  es¬ 
timate  of  rms  float  position  error  is  2.3  meters.  Simulation  results  have  indicated  that  the  true  rms 
error  in  the  least-squares  filter’s  position  estimate  is  no  more  than  about  twice  the  estii;  .ated  error. 
Accordingly,  the  rms  error  in  the  least-squares  filter  estimate  is  thought  to  be  less  than  4.6  meters. 

The  rms  residual  for  the  least-squares  estimates  and  the  rms  innovation  for  the  Kalman  filter 
estimates  are  shown  as  a  function  of  measurement  (Swallow  float  record)  number  in  Figure  3.16. 
Both  rms  residual  and  innovation  contain  large  spikes  around  records  1600  and  1700.  These 
spikes  are  due  to  missing  surface  echo  travel  time  measurements  in  floats  6  and  8’s  data,  and  the 
values  determined  by  the  adaptive  linear  predictor  must  still  contain  large  errors. 

Figure  3.17  shows  the  least-squares  filter’s  float  horizontal  position  estimates  between  records 
1003  and  2120.  The  position  estimates  indicate  that  the  freely  drifting  floats  dispersed  away  from 
the  center  of  the  float  uiangle  with  floats  0  and  I  moving  to  the  northwest,  floats  2  and  3  to  the 
west,  float  4  to  the  southwest,  and  floats  5, 6,  7,  and  8  to  the  southeast.  The  drifting  pattern  was 
probably  due  to  the  complex  water  movement  (including  eddies,  wind-driven  surface  currents, 
and  the  California  Current)  near  the  experiment  site.  The  least-squares  filter’s  float  depth  esti¬ 
mates  are  also  plotted  in  Figure  3.18. 

The  distance  between  the  two  filters’  position  estimates  is  relatively  small  for  all  floats  with  an 
rms  difference  less  than  3  meters  for  most  of  the  experiment.  Figures  3.19  show  the  enlarged  ver¬ 
sion  of  the  Kalman  filter’s  float  depth  estimates  as  a  function  of  record  number.  The  floats  oscil¬ 
late  irregularly  and  slowly,  with  periods  of  30  to  90  minutes,  which  correspond  roughly  to  internal 
wave  periods.  The  slowly  increasing  depth  trend,  except  for  float  0,  is  probably  caused  by  gradual 
compression  of  the  float  as  it  gets  colder. 

The  float  speed  estimates  derived  from  the  Kalman  filter’s  X,  Y,  and  Z  velocity  component  es¬ 
timates  are  plotted  in  Figure  3.20.  The  float  speeds  appear  to  be  rather  unsteady  (i.e.,  high  process 
noise)  throughout  the  experiment.  The  large-scale  speed  oscillation  with  periods  of  4  to  1 2  hours. 
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Figure  3.15  Mean  innovation  and  RMS  error  of  the  Kalman  and  least-squares  filters, 

July  1989  experiment. 
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Figure  3.17  Float  horizontal  displacement  estimates  using  least-squares  filter 
during  the  July  1989  experiment.  The  circles  mark  the  starting  positions. 
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Figure  3.18  Float  depth  estimates  using  least-squares  filter,  July  1989  experiment. 
The  horizontal  axis  of  the  figure  gives  the  time  in  units  of  Swallow  float  record 
number.  Eighty  records  represent  one  hour. 
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Figure  3.19  Float  depth  estimates  using  Kalman  filter,  July  1989  experiment.  The 
horizontal  axes  of  the  figure  give  the  time  in  units  of  Swallow  float  record  number. 

Eighty  records  represent  one  hour. 
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Figure  3.20  Float  speed  estimates  derived  from  the  Kalman  filter’s  X.  Y.  and  Z 
velocity  estimates,  July  1989  experiment.  The  horizontal  axes  of  the  figure  give 
the  time  in  units  of  Swallow  float  record  number.  Eighty  records  represent  one 

hour. 
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depending  on  float  equilibrium  depth,  is  thought  to  be  caused  by  the  tidal  currents.  Average  float 
speeds  for  the  freely  drifting  floats  are  listed  in  Table  3.12. 

3.6  Summary 

In  this  chapter,  two  localization  methods,  generalized  least-squares  filter  and  Kalman  filter, 
were  reviewed  from  the  performance-criteria  point  of  view.  The  inputs  to  both  filters  using  data 
from  the  July  1989  experiment  were  described  in  detail.  The  procedures  for  tuning  the  filters  w-ere 
described,  and  the  results  from  both  filters  were  discussed  and  compared.  The  rms  position  error 
was  estimated  to  be  less  than  3.3  m  for  the  Kalman  filter  and  less  than  4.6  m  for  the  least-squares 
filter.  Due  to  the  high  process  noise  experienced  by  the  floats  during  the  experiment  and  the  quasi¬ 
vertical-line -array  deployment  geometry  [Culver  and  Hodgkiss,  1988],  the  two  filters  performed 
comparably.  The  two  position  estimates  had  an  rms  difference  of  less  than  3  meters  for  most  of 
the  experiment.  Both  methods  appeared  to  be  capable  of  estimating  float  positions  to  within  the 
desired  accuracy  of  one-tenth  of  a  wavelength  at  the  highest  frequency  of  interest  20  Hz  (7.5  m)  in 
order  to  effectively  beamform  the  VLF  acoustic  data  [Steinberg,  1976]. 


Table  3.12  Average  float  speed,  July  1989  experiment. 


Float  Number 

Average 

Speed 

(meters/hour) 

0 

350 

1 

170 

2 

220 

3 

180 

4 

80 

5 

40 

6 

110 

7 

90 

8 

120 
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4.  SOURCE  LOCALIZATION  WITH  SWALLOW  FLOAT  ARRAY 


4.1  Introduction 

This  chapter  presents  the  results  from  matched-field  processing  of  a  14-Hz-tone  propagation 
collected  during  the  July  1989  Swallow  float  experiment  in  the  North-East  Pacific.  There  are  three 
steps  involved  in  the  matched-field  processing  as  depicted  in  Figure  4.1.  The  first  step  is  the  esti¬ 
mation  of  array  covariance  matrix  at  the  source  frequency;  the  second  is  the  prediction  of  replica 
vectors  for  all  assumed  source  locations;  and  the  last  step  is  the  computation  of  ambiguity  surfac¬ 
es:  a  peak  in  the  ambiguity  surfaces  indicates  a  likely  source  location.  The  fundamental  and  theo¬ 
retical  backgrounds  of  MFP  technique  were  given  in  chapter  1. 


Figure  4.1  Matched-field  processing  block  diagram. 


This  chapter  is  organized  as  follows.  Section  2  describes  the  preparation  of  VLF  acoustic  data 
collected  during  the  1989  experiment  and  the  estimation  of  the  array  covariance  matrix  for  the  se¬ 
lected  time  interval.  Section  3  describes  the  acoustic  environment  of  the  oceanic  channel  between 
the  MARK  VI  towed  source  and  the  Swallow  float  array  and  presents  the  modeling  the  replica 
pressure  field.  Section  4  presents  the  MFP  results  on  experimental  data  using  Bartlett.  Minimum 
Variance,  and  MUSIC  processors,  while  section  5  gives  the  controlled  MFP  simulation  to  aid  in 
interpreting  the  MFP  results  from  experimental  data.  Lastly,  section  6  provides  a  summary  of  the 
chapter. 
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4.2  VLF  Acoustic  Data  Preparation 


Although  each  Swallow  float  collects  three  channels  of  geophonc  data  and  one  channel  of  hy¬ 
drophone  data,  only  the  omnidirectional  hydrophone  data  are  studied  and  analyzed  in  this  study. 
The  goal  of  preparing  the  acoustic  data  is  to  process  the  acoustic  time  scries  to  a  form  from  which 
the  cross-spectral  density  matrix  or  array  covariance  matrix  can  be  estimated.  This  section  dis¬ 
cusses  the  essential  steps  needed  to  estimate  the  array  covariance  matrix. 

4,2.1  Aligning  Float  Time  Bases 

Since  each  float  contains  its  own  clock  for  timing,  the  first  step  in  pnxessing  the  acoustic  data 
is  to  align  the  time  bases.  If  we  denote  time  according  to  float / s  clock  as  Tj  and  true  elapsed  time 
as  t.  Then, 

Tj(t)  =  hf  +  Xjt  +  aj  (4.1) 

where  hj  is  the  float  clock  acceleration,  g \j  is  the  float  clock  rate,  and  a,  is  the  float  clock  offset,  all 
of  which  arc  constants  unique  to  each  float.  Note  that  ^  is  dimensionless,  and  the  units  of  hj  and  Uj 
arc  time'1  and  time,  respectively.  Ideally,  ht  and  at  arc  equal  to  0,  and  x,  is  unity.  The  direct  ap¬ 
proach  to  the  time  base  problem  is  to  solve  equation  (4.1 )  for  true  time  t  in  terms  of  Ty  hy  gy,  and 
dj  provided  hy  g and  a^arc  known.  The  float  time  base  can  then  be  aligned  to  true  time.  The  hr 
g  j,  and  Uj  can  be  determined  by  ( 1 )  using  a  satellite  synchronized  clock  on  board  ship  to  frequent¬ 
ly  communicate  with  the  floats  in  real  time  to  get  7ys  and  (2)  fitting  the  Tj’s  with  a  second-order 
curve,  the  coefficients  of  the  fitted  curve  give  the  hy  g;,  and  ay 

An  alternative  approach  to  the  time-base  alignment  problem  is  to  estimate  the  relative  clock 
acceleration,  clock  rate,  and  the  offset,  (/i,  -  hj),  (Xi  -  XjX  and  { a ,  -  ap,  respectively.  Once  these 
quantities  arc  known,  we  may  ch<x)sc  one  float  as  a  reference  whose  time  base  will  not  be  shifted 
and  shift  the  ime  base  of  the  rest  of  the  floats’  time  scries  to  align  them  with  the  reference  float’s 
time  base.  It  has  been  shown  (Culver  ct.  al.,  1989)  that  the  reciprocal  path  travel  times  can  be 
combined  to  estimate  the  relative  clock  acceleration,  clock  rate,  and  the  offset  differences,  (h,  - 
hj),  (x,  -  Xj),  and  (a,  -  ap.  That  is, 

A 7^  ~  (45  x  «) 2  (/?,  -  hp  +  (45  x«)  U,-*,)  +  (a-ap  (4.2) 

where  n  is  the  record  number  (each  record  lasts  45  seconds),  and  A77y  is  approximately  equal  to 
one-half  the  difference  between  the  two  reciprocal  measurements  (  t  (n)  -  t  (//) )/  2 .  As  in  the 
first  approach,  the  curvature,  slope,  and  intercept  of  the  second-order  curve  fitted  to  the  travel  time 
differences  provide  estimates  of  the  differences  between  the  float  clock  accelerations,  rates,  and 
offsets.  In  this  study,  we  have  adopted  the  latter  approach  to  align  the  float  time  bases.  Figure  4.2 


58 


Time  (msec) 


600  900  1200  1500  1800  2100  2400 
Record  Number 
(a) 


600  900  1200  1J00  1800  2100  2400 
Record  Number 
(b) 


Record  Number  Record  Number 

(c)  (d) 


Figure  4.2  (a)  Float  2  listening  to  float  8,  (b)  Float  8  listening  to  float  2,  (c)  Travel  time 
difference  between  floats  2  and  8,  and  (d)  Second-order  fit  of  (c). 


illustrates  the  relative  clock  acceleration,  clock  rate,  and  offset  estimates  between  floats  2  and  8. 
Appendix  C  summarizes  the  coefficients  of  the  second-order  curves  fitted  to  the  travel  time  differ¬ 
ences  for  all  float  pairs. 


4.2.2  Selecting  Data  Records 

To  ensure  the  quality  of  the  array  covariance  matrix  estimate,  the  data  records  need  to  be  qual¬ 
ified  and  selected  with  care.  Two  criteria  used  in  selecting  the  data  record  are  described. 

4.2.2.1  Signal-to-Noise  Ratio 

The  first  criterion  in  selecting  data  records  is  to  observe  pressure  power  spectral  estimates  of 
the  time  series  for  signal-to-noise  ratio  at  the  frequency  of  interest.  The  power  spectra,  computed 
between  0  and  25  Hz  and  plotted  in  Figure  4.3  are  obtained  by  processing  3  minutes  of  data 
(records  1143  to  1146).  The  power  spectra  are  derived  from  incoherently  averaging  28, 50  percent 
overlapped,  512-point  FFTs  (97-mHz  bin  width).  A  Kaiser- Bessel  window  with  a  parameter  of 
2.5,  yielding  a  sidelobe  level  of  -57  dB  [Harris;  1978],  weights  the  data  prior  to  each  FFT.  Power 
values  are  calibrated  in  dB  re  1  jxPa.  The  90  percent  confidence  level  in  these  spectra  is  about  ±1 
dB.  The  14  Hz  was  being  transmitted  during  the  time  when  data  records  1120  through  1160  were 
collected.  A  line  at  14  Hz,  about  10  to  15  dB  above  the  background  noise,  can  be  clearly  seen  in 
all  freely  drifting  floats'  pressure  spectra.  The  high  S/N  at  14  Hz  observed  in  the  power  spectra  il¬ 
lustrates  the  good  quality  of  the  1989  Swallow  float  data  sets. 


4.2.2.2  Spatial  Coherence 

The  second  criterion  in  selecting  the  data  is  to  estimate  the  spatial  coherence  between  floats. 
Coherence  is  a  measure  of  the  causality  and  similarity  between  a  pair  of  signals  received  by  two 
different  sensors.  Given  two  stationary  zero-mean  random  processes  x (t)  and  y(t),  the  coherence, 
or  more  precisely,  the  magnitude-squared  coherence  (MSC)  is  defined  as 


lr,,w|2 


Sxx(f)Syy{f) 


(4.3) 


where  Sxy  (J)  is  the  cross  spectral  density  at  frequency /between  x(t)  and  y(t)  with  power  spectra 
Sr ,  if)  and  Svu  (J)  .  The  MSC  function  evaluated  at /is  a  real  value,  conveniently  normalized  to 
lie  between  zero  and  unity.  High  coherence  of  a  line  frequency  harmonic  among  all  float  pairs  not 
only  indicates  the  signals  originate  from  the  same  source  but  assures  high  array  gain  when  the  in¬ 
dividual  sensor  outputs  are  combined  to  form  a  beamformer.  Figure  4.4  shows  the  magnitude- 
squared  coherence  and  the  phase  difference  between  floats  0  and  2  during  records  1040  -1240, 
while  appendix  D  lists  the  MSC  value  between  all  elements  of  the  freely  drifting  array  during 
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Figure  4.4  (a)  Spatial  coherence  at  1 4  Hz  and  (b)  Phase  difference  between  floats  0  and  2 

at  14  Hz  during  records  1040-  1240. 


record  1145.  The  MSC  functions  are  calculated  by  averaging  over  40.96  seconds  of  data,  128- 
point  FFT’s  with  50  percent  overlap,  providing  31  averages.  The  high  coherence  during  records 
1120  to  1 160  is  a  confirmation  that  the  signal  originates  from  the  same  source,  and  the  exhibition 
of  the  smooth  measure  of  the  phase  differential  qualifies  the  1989  Swallow  float  data  set  for  beam¬ 
forming  application. 

4.2.3  Estimating  the  Array  Covariance  Matrix 

The  problem  of  estimating  the  array  covariance  matrix  is  required  by  all  three  processors  as 
discussed  in  chapter  2.  The  array  covariance  matrix  R  can  be  estimated  from  the  measurement 
data  and  is  given  by 


R  -  i  £  X,X"  (4.4) 

1=1 

where  X .  is  the  vector  of  Fourier  coefficients  across  the  array  at  the  frequency  of  interest /corre¬ 
sponding  to  the  time  snapshot  Frequency  bins  of  0.4-Hz  width  centered  at  14  Hz  were  extracted 
from  128-point  FFT’s  (2.56  seconds  of  data)  using  a  Kaiser-Bessel  (with  a  -  2.5  for  first  sidelobe 
levels  of  -57  dB)  window  and  50  percent  overlap  for  each  of  the  nine  floats  over  a  40.96-second 
data  record.  This  results  in  a  K  =  31  dyad  products  (X.XW,)  being  averaged  for  the  covariance  ma¬ 
trix  estimate.  In  practice,  the  number  of  averages  required  to  produce  a  reliable  covariance  matrix 
estimate  is  about  twice  to  three  times  the  number  of  sensors  in  die  array.  In  our  case,  31  averages 
will  suffice. 

4.2.4  Inverting  the  Array  Covariance  Matrix 

Since  the  matrix  must  be  inverted  for  the  MV  processor,  it  should  be  well  conditioned  and  of 
full  rank.  It  is  well  known  that  as  many  dyads  as  sensors  must  be  averaged  to  ensure  full  rank.  Al¬ 
though  a  large  number  of  snapshots  (i.e.,  31)  used  in  (4.4)  to  compute  the  array  covariance  matrix 
R  ensures  the  invertibility  of  the  covariance  matrix  in  theory,  the  condition  number  of  the  matrix 
given  by  the  ratio  of  the  smallest  to  the  largest  eigenvalue  is  below  10'6,  and  additional  effort  is 
required  [Tran,  1990]. 

The  two  most  commonly  used  techniques  in  handling  the  ill-conditioned  covariance  matrix 
are  the  eigenstructure  decomposition  method  and  the  white-noise  stabilization  method.  The  eigen- 
structure  method  uses  the  spectral  theorem  to  express  the  array  covariance  matrix  R  and  its  in¬ 
verse  as  summations  of  outer  products  of  the  eigenvectors  (assuming  R  is  positive  definite  and 
Hermitian). 


(4.5) 


R  =  iwf 

I  =  1 


1=1  ' 


(4.6) 


If  R  is  not  full  rank,  i.e.,  the  N  sensor  array  covariance  matrix  R  with  K  <  N  nonzero  eigenvalues, 
then  we  can  express  the  inverse  in  (4.6)  as 


R4  =  £  i  v,vf 

i  =  i  * 

R+  is  also  called  the  pseudo-inverse  [Penrose,  1955] 


(4.7) 


Another  approach  to  this  problem  is  white-noise  stabilization,  which  is  a  procedure  of  stabiliz¬ 
ing  the  covariance  matrix  with  spatial  white  noise  to  allow  for  matrix  inversion.  This  method  is 
also  called  "diagonal  loading”  [Carlson,  1988].  In  essence,  this  method  adds  to  the  main  diagonal 
of  R  the  quantity 


7  M 

where  7  is  a  fraction  of  noise  with  typical  values  between  10-3  and  KT1 ,  and  tr  is  the  trace  oper¬ 
ator.  Since  the  white-noise  stabilization  method  requires  less  computation,  it  is  used  in  this  thesis. 
The  covariance  matrix  R  is  stabilized  by  adding  to  the  main  diagonal  the  quantity  KT1  (trR/M). 
Using  a  stabilization  factor  7  of  0.1  corresponds  to  introducing  in  the  system  an  uncorrelated  sen¬ 
sor  noise  10  dB  below  the  average  sensor  power  level  [Tran  and  Hodgkiss,  1991]. 

4.3  Acoustic  Propagation  Modeling 

The  second  step  toward  matched-field  processing  is  the  calculation  of  the  replica  vectors.  To 
predict  the  acoustic  pressure  field  received  by  a  sensor  due  to  an  assumed  source,  one  must  model 
the  acoustic  environment  between  the  source  and  the  sensor.  This  section  describes  the  environ¬ 
mental  data  collected  during  the  1989  experiment;  presents  the  modeling  of  the  acoustic  environ¬ 
ment  between  the  MARK  VI  source  and  Swallow  float  array;  and  provides  the  assumptions  made 
in  computing  the  replica  vectors. 
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4.3.1  Environmental  Data 


As  discussed  in  chapter  2,  a  large  quantity  of  environmental  data  was  taken  during  the  1989 
experiment,  including  CTD  (conductivity,  temperature,  and  depth)  casts,  XBTs  (ship  deployed  ex¬ 
pendable  temperature  profilers),  and  AXBTs  (air  deployed  expendable  temperature  probes).  Plot¬ 
ted  in  Figure  4.5  is  a  series  of  sound  speed  profiles  derived  from  the  CTDs,  the  XBTs  and  the 
AXBTs  taken  along  the  approximate  path  between  the  source  and  the  Swallow  float  array  within 
48  hours  of  the  Swallow  float  experiment.  The  sound  speed  profiles  taken  near  the  source  and 
Swallow  float  deployment  site  are  enlarged  and  plotted  in  Figure  4.6.  The  profiles  are  typical  for 
middle  latitudes  in  the  North  Pacific  Ocean  and  are  characterized  by  [Burdic,  1991]  (1)  a  surface 
layer  that  extends  from  the  surface  to  perhaps  100  m,  (2)  a  layer  that  extends  down  to  perhaps  600 
to  800  m  and  is  characterized  by  a  negative  gradient,  and  (3)  a  deep  layer  where  sound  speed  in¬ 
creases  gradually  with  depth  with  positive  gradient.  All  of  the  profiles  along  the  propagation  path 
exhibit  a  depth  excess  that  is  a  necessary  condition  for  long  range-propagation  [Urick,  1983]. 

The  meso-scale  change  in  the  ocean  temperature  structure  across  range  on  the  order  of  several 
thousand  kilometers  results  in  range-dependent  variations  in  the  sound  speed  profile  as  shown  in 
Figure  4.5.  The  slow  varying  sound  speed  structure  variation,  especially  in  the  upper  water  col¬ 
umn  between  120°  W  and  135°  W,  is  believed  to  be  influenced  by  the  cold  waters  of  the  Califor¬ 
nia  Current  coming  down  from  the  north.  The  sound  speed  structure  between  140°  W  and  150°  W 
remains  relatively  stable;  this  part  of  the  ocean  is  known  to  be  environmentally  benign. 

4.3.2  Modeling  with  Adiabatic  Normal  Modes 

The  adiabatic  approximation  of  the  normal  mode  model  provides  a  simple  solution  to  the 
problem  of  acoustic  propagation  in  a  slowly  range-varying  ocean  (refer  to  section  1. 2.2.1).  In  this 
thesis,  we  use  the  adiabatic  mode  theory  to  model  the  acoustic  pressure  field.  Recall  that  the  adia¬ 
batic  mode  theory  solution  in  a  2-D  environment  is: 


M  exp[i\lm{s)ds 

P(r,z)  =  A  £  um  (zo;0)  um  ( z\r ) - -|== == -  (4.8) 

The  modal  sum  is  over  the  number  of  propagating  modes,  M,  that  exist  between  source  and 
the  receiver.  The  sum  involves  the  depth  eigenfunctions  at  the  source  um  (zo;0) ,  and  the  receiver 
um  ( z:r ) ,  as  well  as  the  horizontal  wave  numbers  function  ( s )  reflecting  the  range  depen¬ 
dence  of  the  medium  along  the  path  between  the  source  and  the  receiver.  In  this  study,  the  imple¬ 
mentation  strategy  for  J (s)  ds  is  to  divide  the  full  range  into  a  number  of  segments;  the 
horizontal  wave  numbers  are  calculated  at  a  discrete  sets  of  ranges  where  sound  speed  measure¬ 
ments  are  available. 
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To  determine  the  number  of  propagating  modes,  M,  for  the  oceanic  waveguide  formed  by  the 
source  and  the  Swallow  float  array,  we  need  to  examine  the  modal  depth  eigenfunctions. 

4.3.2. 1  Depth  Eigenfunctions 

The  mode  depth  functions  um  are  determined  by  the  source  frequency,  the  nature  of  the  sound 
speed  profile,  the  water  depth,  and  the  boundary  conditions  at  the  surface  and  the  bottom.  The 
modal  depth  eigenfunctions  at  a  frequency  of  14  Hz  for  the  environment  near  the  Mark  VI  source 
and  the  Swallow  float  deployment  site  are  computed  with  the  ATLAS  normal  mode  model  [Gor¬ 
don  and  Bucker,  1984].  The  first  30  modal  depth  eigenfunctions  are  gray-leveled  and  displayed  in 
Figure  4.7.  Examining  the  upper  portion  of  the  figure  that  corresponds  to  the  environment  where 
the  source  is  located,  we  see  that  there  are  about  22  modes  trapped  in  the  water  column  and  they 
are  non-bottom  interacting;  also,  for  the  source  depth  at  90  m,  the  first  6  modes  are  very  weakly 
excited.  For  the  low-order  modes  to  be  strongly  excited,  the  source  would  have  to  be  placed  be¬ 
tween  its  turning  points  where  the  mode  peaks  up.  The  lower  portion  of  the  figure  corresponds  to 
the  environment  at  the  array  site;  given  the  water  depth  at  this  location,  non-bottom  interacting 
modes  are  limited  to  the  first  16  modes.  Now  let  us  examine  the  pressure  field  modeled  by  the  adi¬ 
abatic  mode  theory,  i.e.,  equation  (4.8).  Note  that  the  depth  functions  enter  into  the  pressure  field 
calculation  as  product  um  {z0; 0)  um  ( z;r )  where  the  one  depth  function  is  evaluated  at  source 
depth  and  the  other  depth  function  is  evaluated  at  the  receiver  depth.  As  result  of  the  product 
um  um  (z’r) » and  assuming  the  bottom  interacting  modes  are  unable  to  propagate  out  to 
long  range  due  to  bottom  attenuation,  only  the  first  16  trapped  modes  will  be  used  in  the  calcula¬ 
tion  of  the  acoustic  pressure  field. 

To  get  a  visualization  of  the  propagation  of  the  trapped  modes,  we  apply  the  sound  speed  pro¬ 
files  collected  near  the  source  to  a  ray  acoustic  model  RASP  [Palmer  and  Powell,  1989].  Figure 
4.8  shows  a  ray  trace  for  a  source  depth  of  90  meters.  Only  those  rays  that  are  refracted  in  the  wa¬ 
ter  column  are  shown.  The  rays  that  reflected  off  the  bottom  (the  bottom  bounce  rays  that  are 
equivalent  to  the  bottom  interacting  modes)  are  omitted  from  the  ray  trace  since  they  will  not  be 
able  to  propagate  out  to  the  Swallow  float  array  about  2500  km  away  from  the  source  due  to  mul¬ 
tiple  bottom  encounters.  Because  of  the  shallow  source  and  the  depth  excess  in  the  sound  speed 
profile,  the  rays  are  either  refracted-refracted  (RR)  rays  that  are  trapped  in  a  deep  sound  channel 
or  refracted,  surface-reflected  (RSR)  rays  that  do  not  intersect  the  bottom  and  periodically  return 
to  the  sea  surface  on  the  far  side  of  the  convergence  zone  (approximately  50  km),  and  serve  to  ex¬ 
tend  the  ranges  of  good  transmission.  Table  4.1  lists  the  equivalent  source  angles  for  modes 
trapped  in  the  water  column  at  a  frequency  of  14  Hz, 


67 


1480  1500  1520  1540 

Sound  Speed  (m/sec) 


0 


5 


10  15  2< 

Mode  Number 


Sound  Speed  (m/sec) 


10  15  20 

Mode  Number 


Linear  scale 


Figure  4.7  Modal  eigenfunctions  at  frequency  of  14  Hz. 


5000 


43.2.2  Horizontal  Wave  Numbers 


Like  the  mode  depth  eigenfunctions,  the  eigenvalues  or  modal  horizontal  wave  numbers  | 
are  also  determined  by  the  frequency,  the  nature  of  the  sound  speed  structure  (and  hence  water 
depth),  and  the  boundary  conditions  at  the  surface  and  bottom.  In  this  study,  we  use  ATLAS  to 
calculate  the  horizontal  wave  numbers  for  the  discrete  segments  (each  approximately  30  km) 
along  the  path  of  signal  propagation.  Table  4.2  gives  the  eigenvalues  for  three  representative  loca¬ 
tions  along  the  path  between  the  source  and  the  array  for  the  first  16  modes. 


Table  4.1  Equivalent  source  angles  for  modes  trapped  in  the  water  column  at  a  frequency  of  1 4 

Hz. 


Mode  Number 

Source  Angle 
(degree) 

1 

2.0 

2 

5.0 

3 

6.4 

4 

7.5 

5 

8.5 

6 

9.3 

7 

10.0 

8 

10.7 

9 

11.3 

10 

11.9 

11 

12.5 

12 

13.0 

13 

13.5 

14 

13.9 

15 

14.4 

16 

14.8 
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Table  4.2  Eigenvalues  for  modes  trapped  in  the  water  column  at  a  frequency  of  i  4  Hz. 


Mode 

Number 

Eigenvalues 
at  source 
(rad/m) 

Eigenvalues 
at  140°  W 
(rad/m) 

Eigenvalues 
at  SF 
(rad/m) 

t 

0.0593778 

0.0593764 

0.0593557 

2 

0.0591910 

0.0592177 

0.0592227 

3 

0.0590471 

0.0590707 

0.0590927 

4 

0.0589048 

0.0589325 

0.0589622 

5 

0.0587666 

0.0587980 

0.0588299 

6 

0.0586510 

0.0586674 

0.0586971 

7 

0.0585233 

0.0585401 

0.0585663 

8 

0.0583968 

0.0584139 

0.0584366 

9 

0.0582721 

0.0582893 

0.0583082 

10 

0.0581492 

0.0581663 

0.0581831 

11 

0.0580283 

0.0580445 

0.0580601 

12 

0.0579092 

0.0579245 

0.0579396 

13 

0.0577918 

0.0578062 

0.0578219 

14 

0.0576766 

0.0576898 

0.0577127 

15 

0.0575628 

0.0575750 

0.0576024 

16 

0.0574503 

0.0574618 

0.0574595 

4.3.2.3  3-D  Field  Approximation 

As  pointed  out  in  chapter  1 ,  range-dependent  propagation  modeling  using  the  adiabatic  nor¬ 
mal  mode  model  is  a  2-D  implementation  from  which  the  field  at  a  particular  range  and  depth  can 
be  computed,  giving  the  bathymetry  and  the  sound  speed  measurements  along  the  path  between 
the  source  and  receiver.  In  this  thesis,  the  3-D  replica  pressure  field  was  approximated  by  evaluat¬ 
ing  the  adiabatic  model  along  a  series  of  N  bearings  in  the  range-dependent  environment  assum¬ 
ing  the  environment  is  azimuth-invariant,  to  give  an  N  x  2-D  description  of  the  field.  This 
assumption,  which  is  important  for  computational  efficiency  and  assumes  the  sound  paths  in  the 
range-azimuth  plane  are  taken  to  be  linear,  can  nonetheless  provide  a  good  approximation  of  the 
true  3-D  field  [Perkins  and  Baer,  1982]. 
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4.4  Experimental  Data-Processing  Results 


After  aligning  the  float  time  bases  and  performing  data  quality  checks,  the  array  covariance 
matrix  for  record  1145  was  estimated,  and  the  replica  vectors  for  the  matched-  field  processors 
were  generated  using  (4.8).  The  problem  of  grating  lobes  in  beamforming  using  a  sparse  array 
was  not  addressed  in  this  study,  we  therefore  limited  our  focus  to  a  region  of  interest.  The  replica 
vectors  were  thus  computed  for  a  hypothetical  source  in  a  spatial  window  extending  in  range  from 
2300  km  to  2700  km,  in  depth  from  0  to  300  m,  and  in  azimuth  from  166°  to  176°  (refer  to  Fig¬ 
ure  3.14,  we  use  the  mathematical  convention  with  -X  pointing  0°  as  reference  and  rotating  coun¬ 
terclockwise).  The  sampling  intervals  in  range,  depth,  and  azimuth  were  1000  m.  10  m,  and  0.1  °, 
respectively.  Matched-field  processing  was  performed  with  all  three  processors:  Bartlett.  MV,  and 
MUSIC.  The  peak  value  in  the  region  of  interest  was  recorded  and  normalized  to  yield  power  in 
dB  re  pPa.  Table  4.3  summarizes  the  results  of  matched-field  processing  on  the  experimental  data. 


Table  4.3  Matched  field  processing  results. 


Bartlett 

Processor 

MV 

Processor 

MUSIC 

Processor 

Depth  of  Max. 

10m 

130  m 

10m 

Range  of  Max. 

2659  km 

2659  km 

2659  km 

Azimuth  of  Max. 

172.1  ° 

172.1  0 

172.1  0 

Power  of  Max. 

82.7  dB 

73.2  dB 

- 

Note  that  power  levels  were  reported  for  the  Bartlett  and  MV  processors  only  since  the  MUSIC 
method  did  not  yield  the  true  power.  Figures  4.9. 4.10.  and  4.11  present  the  matched-field  ambigu¬ 
ity  surfaces.  The  upper  panel  in  each  figure  was  the  range-azimuth  ambiguity  surface  evaluated  at 
the  depth  where  the  highest  peak  occurred  while  the  lower  panel  was  the  range-depth  ambiguity 
surface  evaluated  at  the  azimuth  where  the  highest  peak  occurred.  In  each  case,  the  surface  was 
normalized  to  its  highest  peak  and  was  marked  with  this  symbol,  *.  For  comparison,  the  true 
source  location  was  marked  with  a  A.  While  mismatch  existed,  all  three  processors  were  in  good 
agreement  except  for  their  depth  estimates.  In  fact,  all  three  lacked  the  depth  resolving  power.  The 
highest  peak  in  the  ambiguity  surface  differs  from  the  true  location  by  0.6°  in  azimuth  and  166  km 
in  range.  The  large  number  of  sidelobes  observed  in  the  ambiguity  surfaces  was  thought  to  be  due 
to  imperfect  modeling  and  the  sparseness  of  the  array.  The  MV  and  MUSIC  processors  suffered 
large  losses  due  to  mismatch  since  the  replica  were  imperfect.  Ambiguities  in  range,  which  have 
been  referred  to  as  sidelobes,  were  the  result  of  the  repetitive  structure  of  the  acoustic  field  in  a 
convergence  zone  environment. 


Azimuth  (degree) 


Depth  (m)  Azimuth  (degree) 


Range  (km) 


-3  -2  -I  0 

Relative  Power  (dB) 


Figure  4.10  Matched-field  processing  results  using  the  MV  method  for  the 
14-Hz  source  during  record  1145. 
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Figure  4.1 l  Matched-field  processing  results  using  the  MUSIC  method  for 
the  1 4-Hz  source  during  record  1145. 
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In  addition  to  the  localization  performance  computed  in  Table  4.3,  the  array  signal  gain  was 
also  investigated.  Several  array  performance  measures  for  matched-field  processing  have  been 
proposed  in  the  literature.  In  this  thesis,  we  define  array  signal  gain  as  the  ratio  of  matched-field 
processor  output  S/N  divided  by  the  average  input  S/N  on  the  floats.  However,  the  signal  level 
was  expected  to  be  nonuniform  over  the  array,  so  an  alternative  choice  of  denominators  such  as 
median  float  input  S/N  (signal  level)  and  maximum  float  input  S/N  (signal  level)  were  also  com¬ 
puted.  Table  4.4  lists  the  array  gain  for  various  choices  of  definition  using  the  Bartlett  processor 
for  record  1145.  Note  that  the  theoretical  array  gain  assuming  white  noise  at  the  sensor  output  for 
a  nine  sensor  array  is  10  log9  =  9.54. 

Table  4.4  Matched-field  processing  array  gain. 


Bartlett  Processor 

Array  Gain 
(dB) 

Max.  input  S/N 

2.6 

Ave.  input  S/N 

5.9 

Med.  input  S/N 

6.8 

4.5  Controlled  Simulation 

While  source  localization  in  azimuth  was  somewhat  successful,  localization  in  range  and 
depth  seemed  to  be  a  problem.  Simulations  are  presented  here  in  an  effort  to  understand  the  exper¬ 
imental  ambiguity  surfaces.  Three  simulation  cases  were  studied:  the  ideal  simulation,  uncertainty 
in  float  positions,  and  uncertainty  in  sound  speed  structure. 

4.5.1  No  Mismatch  Simulations 

Assuming  there  was  no  mismatch  and  input  S/N  was  10  dB,  the  simulated  "acoustic  data”  and 
replica  vectors  were  generated  using  the  same  environment  model.  A  14-Hz  source  was  simulated 
at  a  range  of  2493  km,  an  azimuth  of  171.5°,  and  a  depth  of  90  m,  which  is  the  true  source  loca¬ 
tion  at  record  1145  according  to  the  R/V  Aloha  source  log.  The  ambiguity  functions  for  all  3  pro¬ 
cessors  were  evaluated  at  a  depth  of  90  m  and  an  azimuth  of  171.5°  degrees  and  plotted  in  Figures 
4.12,  4.13,  and  4.14,  respectively.  As  expected,  the  source  was  correctly  localized.  Since  there 
was  no  mismatch,  the  dynamic  range  of  the  ambiguity  surfaces  was  much  larger  than  in  the  case 
of  the  real  data.  The  Bartlett  processor  produced  a  large  number  of  sidelobes;  the  MV  processor 
had  better  suppression  of  sidelobes  and  higher  resolution  of  the  main  lobe:  and  the  MUSIC  pro¬ 
cessors  produced  a  single  peak  in  the  range-depth  cell  of  the  simulated  source. 
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Figure  4.13  Matched-field  processing  simulation  results  (no  mismatch) 

using  the  MV  method. 
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Figure  4.14  Matched-field  processing  simulation  results  (no  mismatch) 

using  the  MUSIC  method. 
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The  high  sidelobes  found  predominately  for  the  ambiguity  surfaces  produced  by  the  Bartlett 
processors  were  believed  to  be  due  to  the  nature  of  the  processor  and  the  array  geometry.  Also,  the 
pressure  field  calculated  using  the  adiabatic  normal  mode  model  was  composed  of  only  the  first 
16  low-order,  water-borne  modes.  Thus,  the  acoustic  pressure  field  was  less  complex  or  less 
unique  at  the  simulated  source  location. 

The  poor  depth  resolution  observed  in  the  Bartlett  and  MV  ambiguity  surfaces  was  thought  to 
be  due  to  the  combination  of  few  propagating  modes  and  the  low  source  depth  (90  m)  to  wave¬ 
length  (105  m)  ratio.  Not  too  surprisingly,  the  MUSIC  processor  yields  excellent  depth  resolution 
due  to  exploitation  of  the  orthogonality  principle. 

4.S.2  Uncertainty  in  Float  Positions 

The  sensitivity  to  sensor  position  mismatch  was  then  investigated.  Assuming  there  were  no 
position  eirors,  the  replica  field  were  generated  using  the  same  conditions  as  in  the  real  data  case. 
However,  the  simulated  “acoustic  data”  was  computed  using  the  array  geometry  that  incorporates 
random  perturbations  of  7  meters  to  each  sensor  position.  The  impact  of  mismatch  due  to  float  po¬ 
sition  error  was  investigated  again  in  a  10-dB  input  S/N  case.  The  ambiguity  surface  is  plotted  in 
Figures  4.15,  4.16  and  4.17  and  the  matched  field  processing  results  are  summarized  in  Table  4.5. 


Table  4.5  Matched-field  simulation  of  sensor  position  errors. 


Bartlett 

Processor 

MV 

Processor 

MUSIC 

Processor 

Depth  of  Max. 

110  m 

110m 

110m 

Range  of  Max. 

2493  km 

2493  km 

2493  km 

Azimuth  of  Max. 

171.7° 

171.7° 

171.7° 

Power  of  Max. 

-0.16  dB 

-7.25  dB 

- 

The  sidelobe  structures  were  very  similar  to  those  of  ideal  simulation,  but  the  peak  value  for 
the  MV  and  for  the  MUSIC  processor  was  much  reduced  from  that  of  the  ideal  simulation.  Al¬ 
though  the  mismatch  reduces  the  dynamic  range  of  the  MV  and  MUSIC  processors  tremendously, 
the  source  was  successfully  located  in  range  and  azimuth  with  reduced  power.  The  estimated 
source  range  was  identical  to  the  “true”  source  range  with  a  minor  discrepancy  in  the  azimuth  es¬ 
timate.  Unfortunately,  the  MUSIC  processor  failed  to  maintain  the  depth  resolving  power  as  in  the 
ideal  case.  These  simulated  results  suggest  that  slight  float  position  error,  i.e.,  less  than  one  tenth 
of  the  wavelength  (10.7  m)  might  be  the  cause  of  the  mismatch  in  depth  and  azimuth  but  cannot 
be  held  responsible  for  the  mismatch  in  range. 
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4.5.3  Uncertainty  in  Sound  Speed  Structure 

We  next  studied  the  simulation  of  sound  speed  mismatch.  The  sound  speed  profiles  collected 
during  the  experiment  (refer  to  Figure  4.5)  were  modified  by  adding  a  linear  function  of  the  form 
[Porter  et.  al.,  1987]: 


Ac(z) 


5  (D-z) 
D 


meter/second 


(4.9) 


to  the  sound  speed  profiles  collected  between  140°  W  and  150°  W.  In  this  simulation,  we  used 
6  =  -3  and  D  -  2000  m  so  that  at  the  surface,  the  sound  speed  was  decreased  by  3  meters/sec¬ 
ond,  while  below  2000  m,  there  was  no  change.  This  modification  was  justified  by  the  fact  that  the 
sound  speed  profiles  collected  in  this  track  were  off  the  signal  propagation  path  (refer  to  Figure 
2.6).  The  replica  vectors  were  generated  using  the  original  profiles  while  the  simulated  "acoustic 
data”  were  generated  using  the  modified  sound  speed  profiles  reflecting  (4.9). 

The  ambiguity  surfaces  for  the  sound  speed  mismatch  simulation  are  plotted  in  Figures  4.18. 
4.19  and  4.20.  and  the  matched-field  processing  results  are  summarized  in  Table  4.6.  The  ambigu- 

Table  4.6  Matched-field  simulation  of  sound  speed  mismatch. 


Bartlett 

Processor 

MV 

Processor 

MUSIC 

Processor 

Depth  of  Max. 

110 

110 

110 

Range  of  Max. 

2610 

2610 

2610 

Azimuth  of  Max. 

171.4° 

171.4° 

171.4° 

Power  of  Max. 

-0.5  dB 

-10.7  dB 

- 

ity  surfaces  show  the  source  peak  shifted  in  range  with  much  reduced  power  for  the  MV  and  MU¬ 
SIC  processors.  The  estimated  source  range  is  off  by  117  km  from  the  “true”  source  range,  and  the 
source  azimuth  is  slightly  shifted  by  0.1°.  These  simulated  results  confirmed  that  uncertainty  in 
sound  speed  structures  can  be  held  responsible  for  the  large  mismatch  observed  in  the  real-data 
ambiguity  surfaces,  particularly  the  range  error. 

4.6  Summary 

This  chapter  presented  the  results  from  the  matched-field  processing  of  data  from  14-Hz  cw 
collected  by  MPL’s  Swallow  float  array  in  a  1989  experiment.  The  data  used  to  estimate  the  array 
covariance  matrix  were  selected  at  time  intervals  during  which  the  signal-to-noise  ratio  at  the  out- 
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Range-Azimuth  ambiguity  surface  at  depth  110  m  (MV) 


Relative  Power  (dB) 


Figure  4.19  Matched-field  processing  simulation  of  uncertainty  in  sound 
speed  sUucture  using  the  MV  methou. 


Figure  4.20  Matched-field  processing  simulation  of  uncertainty  in  sound 
speed  structure  using  die  MUSIC  method. 


put  of  each  individual  float  was  high  and  the  magnitude-  squared  coherence  at  source  frequency 
was  high  among  all  floats.  Eighty-four  sound  speed  profiles  derived  from  AXBT,  XBT,  and  CTD 
measurements  collected  between  34°  50'  N,  122°  20'  W  and  32°  N,  150°  W  within  48  hours  of  the 
Swallow'  float  experiment  were  input  to  an  adiabatic  normal  mode  model  to  compute  the  replica 
vectors.  In  spite  of  the  presence  of  high  sidelobes  and  grating  lobes  due  to  the  sparse  array,  the 
matched-field  ambiguity  surfaces  generated  by  all  three  MFP  processor  were  consistent.  While 
the  azimuth  estimate  w>as  quite  successful,  all  three  processors  experienced  difficulties  «n  resolv¬ 
ing  the  source  depth,  and  all  three  exhibited  errcs  in  range  estimates  by  three  convergence  zones. 
Controlled  simulations  were  performed  to  aid  in  interpreting  the  real  data  processing  results  and 
confirmed  that  (1)  depth  resolution  is  indeed  a  difficult  problem  due  to  few  modes  propagating  out 
to  long  range  and  due  to  a  low  source  depth-to-wavelength  ratio,  (2)  the  range  estimate  is  sensi¬ 
tive  to  uncertainty  in  sound  speed  structure,  and  (3)  the  azimuth  estimate  is  robust.  The  environ¬ 
mental  mismatch  problem  will  be  further  investigated  in  chapter  5. 
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5.  SELF-COHERENT  MATCHED-FIELD  PROCESSING 


5.1  Introduction 

Matched-field  processing  has  been  developed  for  localizing  underwater  acoustic  sources  by 
comparing  acoustic  data  with  predicted  replica  pressure  fields.  The  inputs  to  MFP  are  the  acoustic 
parameters  of  the  ocean  for  predicting  the  replica  pressure  fields  and  the  acoustic  data  received  by 
an  array  of  hydrophones.  In  order  to  compute  the  replica  vectors  for  the  matched-field  processor, 
the  knowledge  of  sound  speed  structure  and  bottom  characteristics  as  a  function  of  depth  and 
range  are  required.  As  shown  in  the  previous  chapter,  if  the  available  environmental  information 
is  not  sufficiently  accurate  (a  situation  referred  to  as  mismatch)  MFP  can  be  degraded  even  if  the 
signal-to-noise  ratio  is  high.  As  an  added  complication,  time-varying  oceanographic  features  such 
as  meso-scale  eddies,  ocean  fronts,  and  internal  gravity  waves  create  fluctuations  that  result  in  the 
time-varying  defocussing  of  the  matched-field  processor.  Thus,  calibration  of  the  environmental 
parameters  so  as  to  improve  the  matched-processing  performance  is  of  special  interest.  This  chap¬ 
ter  investigates  such  an  issue  and  is  organized  as  follows.  Section  2  reviews  two  approaches  to  the 
environmental  mismatch  problem  reported  in  the  literature.  Section  3  proposes  an  approach  to  the 
environmental  mismatch  problem.  Section  4  discusses  the  processing  results  using  the  proposed 
technique.  Section  5  presents  the  chapter  summary. 

5.2  Literature  Review 


5.2.1  Self-Cohering  Technique 

The  self-cohering  technique  proposed  by  Bucker  [Bucker,  19S9;  Tran  and  Hodgkiss,  1992]  in¬ 
volves  using  a  narrowband  low  frequency  source  at  a  known  location  in  an  imprecisely  known 
oceanic  environment.  The  essence  of  the  self-cohering  technique  is  to  introduce  a  correction  fac¬ 
tor  to  the  matched-field  correlation  formulation  to  compensate  for  the  uncertainty  in  the  ocean- 
acoustic  environment  so  that  the  matched-field  correlation  between  the  received  pressured  fields 
and  the  replica  vector  is  maximized.  The  correction  factor  has  shown  to  be  range-dependent  and 
thus  can  be  used  to  correct  for  the  environmental  parameters  (i.e.,  the  modal  phases  and  the  modal 
amplitudes)  in  a  subsequent  matched-field  localization  in  search  for  other  unknown  source(s). 

Recall  that  the  modeled  pressure  field  in  a  weakly  range-dependent  environment  using  adia¬ 
batic  normal  mode  modeling  is  expressed  as 


p  ( r ,  z) 


M 


m  ~  1 


eXP(ij%mdr') 


(5.1) 
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where  z  is  the  depth  of  the  sensor,  zs  is  the  depth  of  the  source,  r  is  the  range  of  the  source,  um  is 
the  m,h  mode  eigenfunction,  is  the  mth  mode  horizontal  wave  number,  and  M  modes  are  used 
to  model  the  acoustic  pressure  field. 

If  we  express  the  Bartlett  matched-field  processor  output  power  as 


P  = 


Y/ih 


(5.2) 


where  J  is  the  total  number  of  sensors,  pj  is  the  complex  pressure  measured  at  the  }h  sensor,  and 
the  pj  is  the  predicted  complex  pressure  at  \he/h  sensor.  Equation  (5.1)  can  be  substituted  into 

(5.2),  and  the  summations  on  the  sensors  and  on  the  modes  can  be  rearranged  to  yield 


P  = 


(5.3) 


where  the  first  term  corresponds  to  the  receiving  part 


and  the  second  term  corresponds  to  the  source  part 


(5.4) 


=  w. 


(v°) 


exp 


Jr 


(5.5) 


Correction  factors  for  each  mode  Cm  are  introduced  into  the  summation  of  Equation  (5.3) 


P  = 


(5.6) 


to  maximize  the  correlation  for  each  mode  with  respect  to  a  source  at  known  location  (ro,  zQ) . 
Using  this  criterion,  the  correction  factor  Cm  is  given  at  the  known  source  location  by 


C 


m 


(5.7) 


where  the  known  source  range  and  depth  (ro,  z0)  are  used  in  Equation  (5.5)  to  compute  Sm.  The 
correction  factors  are  a  function  of  range: 
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(5.8; 


C„M)  =  \Cm(ro'>\exP(i'-ar8(Cm(r0))) 

'  O 

It  has  been  demonstrated  in  simulation  that  (5.8)  calculated  using  a  source  at  a  known  location  can 
be  used  in  search  of  other  source(s)  at  a  longer  range. 

5.2.2  Focalization  Technique 

Focalization  [Collins  and  Kuperman,  1991]  is  an  approach  for  reducing  the  problem  of  mis¬ 
match  by  treating  the  ocean  as  an  acoustic  lens  that  is  out  of  focus  and  by  adjusting  the  acoustic 
parameters  until  the  source  come  into  focus.  This  approach  does  not  require  a  reference  source  but 
assumes  the  measured  modal  phases  are  available,  although  determining  the  measured  modal 
phases  can  be  difficult  as  pointed  out  by  [Collins  and  Kuperman,  1991].  The  focalization  is  per¬ 
formed  by  backpropagating  the  modeled  modal  phases  using  adiabatic  normal  mode  approxima¬ 
tion  until  the  adiabatic  modal  phases  match  the  measured  modal  phases.  The  simulated  annealing 
optimization  technique  is  used  to  search  the  modal  phase  parameter  spaces  via  sound  speed  struc¬ 
ture  for  the  optimal  values.  It  has  been  demonstrated  by  the  focalization  technique  that  the  ocean 
parameters  that  bring  the  source  into  focus  may  not  be  unique;  the  focalization  parameter  search 
sometimes  converges  to  the  correct  source  location  without  converging  to  the  correct  environ¬ 
mental  parameters. 

5.3  Environment  Adaptation  Technique 

In  this  section,  we  present  an  approaih  to  the  environment  mismatch  problem.  We  envision 
that,  similar  to  the  self-cohering  technique  by  Bucker,  the  way  that  the  environmental  mismatch 
can  be  reduced  is  a  two-phase  process: 

1.  Adaptation  phase:  During  this  phase,  a  narrowband  signal  with  frequency  of  interest  at  a 
known  location  is  transmitted  to  probe  the  oceanic  waveguide.  The  signal  could  be  a  sur¬ 
face  ship  of  opportunity  or  a  broadband  source  with  good  S/N  on  the  frequency  of  interest 
such  as  air-deployed  shots.  The  matched-field  processor  is  configured  in  a  feedback  loop 
fashion  as  in  Figure  5.1  to  adjust  the  environmental  parameters  with  the  goal  of  causing 
the  predicted  pressure  field  to  match  the  measured  pressure  field,  i.e., 

maximize  zo’ (^.9) 

r 

where  T  is  the  environmental  parameter  set,  and  P(ra,  z0,  0O)  is  the  matched-field  pro¬ 
cessor  output  power  due  to  source  at  a  known  location  ( ro ,  z0,  BJ  . 
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2.  Localization  phase:  When  the  environment  adaptation  phase  is  completed,  the  optimized 
environmental  parameter  set  ropl  is  then  used  to  compute  the  replica  pressure  fields  and 
normal  matched-field  processing  can  resume  to  search  for  an  unknown  target  of  interest  in 
the  vicinity  of  the  reference  source. 


Figure  5.1  Environment  adaptation  technique  block  diagram. 


In  this  study,  two  types  of  environmental  parameters,  sound  speed  structure  and  horizontal 
wavenumber,  are  considered.  We  use  matched-field  processor  output  as  a  performance  function  or 
cost  function,  since  the  cost  function  P  is  nonlinear  and  may  have  many  local  maxima  (or  local 
minima  if  we  use  the  negative  of  P  as  the  cost  function);  a  global  optimization  technique  is  re¬ 
quired  for  searching  the  optimum  environmental  parameters.  This  section  presents  the  global  op¬ 
timization  method  while  the  processing  results  using  the  environment  adaptation  technique  are 
given  in  the  next  section. 

5.3.1  Global  Optimization  Method 

In  this  section  we  present  an  efficient  global  optimization  method  for  searching  the  environ¬ 
mental  parameter  spaces.  The  method  is  based  on  the  conventional  simulated  annealing  tech- 
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nique.  A  detailed  description  of  the  technique  is  given  in  [Kirkpatrick  et.  al.,  1983].  For 
completeness,  the  motivation  and  a  brief  description  of  the  algorithm  are  given  here. 

5.3. 1.1  Motivation 

Traditionally,  optimization  problems  are  often  treated  by  the  iterative  improvement  method 
such  as  gradient  search.  Iterative  improvement  starts  at  some  initial  state,  or  set  of  values  for  the 
variables.  Then  successive  small  changes  in  the  state  are  made  in  such  a  way  that  the  cost  function 
is  decreased  at  each  step.  The  moves  from  one  state  to  another  are  effected  by  some  procedure. 
This  process  continues  until  there  are  no  moves  that  reduce  the  cost  function.  At  this  point  we 
have  found  a  minimum.  Often,  however,  it  is  a  local  rather  than  a  global  minimum  if  the  perfor¬ 
mance  surface  of  the  functional  is  multi-modal.  A  direct  response  to  this  problem  is  to  perform  it¬ 
erative  improvement  from  several  initial  states  and  select  the  best  of  the  resulting  local  minima. 
The  simulated  annealing  algorithm  is  a  more  effective  rrrthod  that  avoids  this  pitfall  by  allowing 
•‘hill  climbing”  moves  with  probability.  That  is,  it  may  accept  moves  that  generate  a  cost  function 
with  a  higher  cost  than  the  current  one;  thus  it  can  pull  itself  out  of  the  local  minima  and  potential¬ 
ly  fall  into  a  more  promising  down-hill  path. 

The  simulated  annealing  is  suggested  by  an  analogy  between  the  search  for  a  minimum  in  a 
cost  function  and  the  physical  process  by  which  a  material  changes  states  while  minimizing  its  en¬ 
ergy.  In  this  analogy,  a  random  initial  state  (or  set  of  parameter  values)  in  the  optimization  corre¬ 
sponds  to  the  melted  state,  and  local  minima  in  the  cost  function  relate  to  the  formation  of 
different  crystal  structures  once  the  material  freezes  into  a  solid.  In  metallurgy,  a  technique  called 
annealing  is  often  used,  where  a  very  fine  crystal  is  obtained  through  gradual  cooling.  Too  rapid  a 
cooling  schedule  results  in  a  coarse  structure,  similar  to  a  poor  local  minimum  in  an  optimization. 
The  correspondence  between  the  two  processes  was  first  introduced  by  [Metropolis  et.  al.,  1953]. 

5.3.1. 2  Simulated  Annealing  Algorithm 

For  a  system  to  be  optimized,  a  cost  function  £  is  first  established  and  a  dynamic  variable  T, 
the  “temperature”  of  the  system,  is  chosen  to  control  the  process.  Starting  at  a  high  temperature, 
the  system  is  slowly  cooled  down  until  the  system  “freezes”  and  reaches  the  optimum  state  m  a 
manner  similar  to  the  annealing  of  a  crystal  during  growth  to  reach  a  near-perfect  structure.  At 
each  “temperature,”  a  change  in  the  system  state  is  made  according  to  a  certain  rule,  and  then  the 
“energy”  or  “cost”  change  of  the  system  A E  is  calculated.  If  A£  <  0,  the  system  state  alteration  is 
accepted,  and  the  process  is  ctgjLinued.  If  A£  >  0,  then  the  system  state  alteration  is  accepted  with 
probability  defined  as  exp  (~jf)  >  the  well-known  Boltzmann  probability  distribution  in  thermo¬ 
dynamics.  In  the  optimization  problem,  the  quantity  k  is  a  constant,  and  its  dimension  depends  on 
the  dimensions  of  A  £  an^T.  Then  a  random  number  %  uniformly  distributed  in  the  interval  [0,1] 
is  chosen.  If  x  -  exP  (~Jf )  * the  change  of  system  state  is  accepted;  on  the  other  hand,  if 
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X  >  exp  (— — ) ,  the  change  is  discarded,  that  is,  the  system  state  before  change  is  used  for  the 
next  step  of  the  process.  This  procedure  is  repeated  for  each  temperature  until  the  system  is  opti¬ 
mized  by  arriving  at  a  global  energy  minimum.  The  choices  of  initial  temperature  T()  and  the 
cooling  schedule  are  crucial  for  the  success  and  the  speed  of  convergence  of  the  simulated  anneal¬ 
ing  process.  Due  to  the  probabilistic  Boltzmann  selection  rule  of  the  A E  >  0  case,  the  process  can 
always  get  out  of  a  local  minimum  of  the  cost  function  in  which  it  could  get  trapped  and  proceed 
to  the  desired  global  minimum.  This  makes  simulated  annealing  different  from  the  iterative  im¬ 
provement  procedure.  To  make  use  of  the  annealing  algorithm  for  the  optimization  problem,  one 
needs  to  provide  the  following  basic  components  [Press  et.  al.  1988]: 

1.  A  cost  function  E  whose  minimization  is  the  goal  of  the  procedure. 

2.  A  model  of  what  a  legal  set  of  parameter  values  is.  This  model  represents  the  possible 
problem  solutions  over  which  we  will  search  for  a  good  answer. 

3.  Generators  of  random  changes  in  the  system  state  to  determine  how  many  variables  and 
which  variables  will  be  altered  and  how  many  perturbations  there  will  be  to  the  current 
values  of  the  selected  variables. 

4.  A  control  parameter  T  and  an  annealing  (or  cooling)  schedule.  Specifically,  we  need  a 
starting  hot  temperature  (or  a  heuristic  for  determining  a  starting  temperature  for  the  prob¬ 
lem)  and  rules  to  determine  when  the  current  temperature  should  be  lowered,  by  how 
much  the  temperature  should  be  lowered,  and  when  annealing  should  be  terminated. 

5.3. 1.3  Fast  Simulated  Annealing  Procedure 

Our  implementation  of  the  simulated  annealing  algorithm  simply  uses  the  negative  of  the 
matched-field  output  power  as  the  energy  (or  cost  function)  to  be  minimized,  and  sound  speed  or 
modal  wave  number  EOF  coefficients  as  the  environmental  parameters  to  be  varied  (refer  to  sec¬ 
tion  5.3.2)  during  the  annealing  process.  We  use  a  random-number  generator  with  Poisson  distri¬ 
bution  to  determine  the  numbei  of  parameters  to  be  changed;  a  random-number  generator  with 
uniform  distribution  to  select  the  parameters  to  be  varied;  and  a  random-number  generator  with 
Cauchy  distribution  to  generate  the  perturbations  for  the  selected  variables.  The  merit  of  the 
Cauchy  distribuiion  is  that  it  has  a  Gaussian-like  peak  and  Lorentzian  tails  that  imply  more  often 
large  perturbations  in  the  step  space  among  dense  small  perturbations  around  the  current  parame¬ 
ters)  that  allow(s)  fast  escape  from  the  local  minima  and  fast  convergence.  As  a  result  of  using 
the  Cauchy  generating  density,  we  adopted  a  fast  annealing  schedule  suggested  by  [Szu  and  Han¬ 
ley,  1987]  which  is  inversely  linear  in  time,  i.e.,  T-  =  T0f  i,  where  the  T0  is  the  initial  starting 
temperature,  and  i  is  the  iteration  number  (note  that  the  conventional  simulated  annealing  algo¬ 
rithm  uses  a  Gaussian  distribution  to  generate  the  perturbations,  and  it  has  been  proven  by  |Ge- 
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man  and  Geman,  1987)  that  the  cooling  schedule  must  be  inversely  proportional  to  the  logarithm 
of  iteration  /  in  order  to  guarantee  convergence  to  a  global  minimum,  i.e.,  Tt  -  T(>  /  log  i).  7  is 
determined  empirically  to  give  good  results.  In  practice,  the  starting  temperature  is  chosen  to  give 
a  high  accept  rate  when  A£  >  0.  Lastly,  the  stopping  criterion  is  to  terminate  annealing  when  the 
cost  improvement  seen  across  successive  temperatures  is  sufficiently  small.  The  fast  simulated 
annealing  algorithm  for  this  study  is  encapsulated  in  the  following  procedure: 

1.  Set  the  initial  temperature  T0  to  a  large  value. 

2.  Determine  the  number  of  parameters  needed  to  be  perturbed  using  Poisson  distribution 
with  a  selection  rate  of  1. 

3.  Generate  the  perturbations  for  the  selected  parameter(s)  using  Cauchy  distribution. 

4.  Perturb  the  current  parameter  set  and  compute  the  new  cost,  Ej. 

5.  If  £,  <  £(  _  , ,  then  accept  the  new  parameter  set  and  proceed  to  (7). 

6.  If  £(  >£,._,,  then  accept  the  new  parameters  with  a  probability  given  by  the  Boltzmann 
distribution. 

7.  Decrease  the  temperature  according  to  the  fast  cooling  schedule  and  repeat  the  procedure: 
terminate  the  procedure  when  no  new  parameters  are  accepted  for  a  large  number  of  Inter¬ 
actions. 

5.3,2  Reducing  the  Parameter  Search  Spaces 

For  optimization  methods,  an  efficient  parameterization  to  reduce  the  parameter  search  space 
leads  to  fast  convergence  and  more  uniqueness  in  the  solution.  Thus  we  would  like  to  characterize 
the  environment  of  interest  in  as  few  parameters  as  possible  yet  in  a  meaningful  way.  A  common 
method  from  statistics  for  analyzing  data  is  principal  component  analysis.  The  essence  of  this 
method  is  to  find  a  set  of  K  orthogonal  vectors  in  data  space  that  account  for  as  much  as  possible 
of  the  data’s  variance.  Projecting  the  data  from  their  original  A-dimensional  space  onto  the  £-di- 
mensional  subspace  spanned  by  these  vectors  then  performs  a  dimensionality  reduction  that  often 
retains  most  of  the  intrinsic  information  in  the  data.  Topically,  K  «  N,  thus  making  the  reduced 
data  much  easier  to  handle. 

Similar  to  the  principal  component  analysis  method,  oceanographers  have  developed  a  meth¬ 
od  for  deriving  efficient  basis  functions,  known  as  empirical  orthogonal  functions  (EOFs).  for 
measured  physical  quantities  such  as  temperature,  salinity,  or  sound  speed  as  a  function  of  depth 
[Davis,  1976;  Tolstoy  et.  al„  1991).  In  an  effort  to  reduce  the  ocean-acoustic  parameter  spaces. 
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one  can  describe  the  parameter  set  as  a  sum  of  EOFs.  The  EOFs  are  defined  as  the  eigenvectors  V, 
of  the  parameter  covariance  matrix  R: 


RV  =  XV,  (5.10) 

where  X,  is  the  i!h  eigenvalue  or  the  variance  associate  with  V,.  The  covariance  matrix  for  the  en¬ 
vironmental  data  R  is  defined  as 


R  =  E[ffr]  (5.11) 

where  f  =  T  -  E  [F] ,  and  T  is  the  measured  environmental  parameter  values.  E  [  ]  is  the  expec¬ 
tation  (average)  operator. 

Any  one  of  the  environmental  parameter  sets  used  in  estimating  the  covariance  matrix  can  be  rep¬ 
resented  or  spanned  as  a  linear  combination  of  appropriate  eigenvectors: 

N 

f  -  E[n  +  £  a„Vn  (5.12) 

n  =  1 

where  N  is  the  total  number  of  eigenvalues,  and  the  Vn ’s  are  indexed  so  that  Xn  >  Xn  + , .  The  mer¬ 
it  of  using  the  EOFs  to  reduce  the  parameter  search  spaces  is  twofold: 

1.  If  we  express  the  true  parameter  set  in  terms  of  some  functional  expansions: 

r*=E[ritI|i.F,  (513) 

n  =  1 

then  for  a  fixed  number  of  functions  K  <  N,  no  approximate  representation  of  the  parame¬ 
ter  set 


K 

f  =  Etn+XP*F*  (5.14) 

*=  i 

can  produce  a  lower  mean -square  error  e[  (T  -  f)  r(rinje  -  f )]  than  is  obtained  us¬ 
ing  F*  =  \k  according  to  the  least-squares  principle  [Davis,  1976].  In  practice,  a  high  de¬ 
gree  of  accuracy  can  be  achieved  with  only  two  or  three  EOFs  for  ocean -acoustic 
parameters  [LeBlanc  and  Middleton,  1980]. 
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2. 


The  covariance  matrix  calculated  from  measurements  is  a  complete  summary  of  the  man¬ 
ner  in  which  the  environment  varies.  This  matrix  reveals  the  inter-statistical  dependence 
between  parameters.  Thus  the  spanned  parameter  set  is  meaningful  since  the  physics  gov¬ 
erning  the  process  that  generates  the  parameter  set  is  preserved. 

Using  the  EOF  approach  to  parameterize  the  environment,  (5.9)  can  now  be  expressed  as 

maximize  P  (r<?>  zo'  (5.15) 

{ak,  k-l,K} 

where  is  related  to  T  by 

K 

T  =  £[T]+  5>*v*  (5.16) 

*=  i 


5.4  Processing  Results 

The  environment  adaptation  technique  described  above  is  first  applied  to  simulation  data  and 
then  used  to  process  data  collected  with  the  Swallow  float  array  that  was  deployed  in  the  North¬ 
east  Pacific  in  July  1989. 

5.4.1  Simulation 

To  make  the  simulation  as  realistic  as  possible,  we  modeled  the  source-array  geometry  corre¬ 
sponding  to  the  July  1989  experiment  with  a  14-Hz  source  deployed  at  the  locations  according  to 
the  R/V  Aloha  source  log.  The  freely  drifting  sensor  array  geometries  corresponding  to  the  navi¬ 
gation  results  during  the  three  time  intervals,  i.e.,  01:47,  04:38,  and  07:32  PST  July  9,  1989  were 
used  in  the  simulation  study. 

To  simulate  the  environmental  mismatch,  similar  to  that  of  the  last  chapter,  we  modified  the 
profiles  between  34°  N,  140°  W  and  32°  N,  150°  W  by  adding  a  function  of  the  form: 

AC(z)  =  C34.A,>l40.w(z)  “Cmean(z)  meter/second  (5.17) 

where  C34<w  l40„w  is  the  sound  speed  profile  taken  at  34°  N,  1 40°  W,  and  Cmean  is  the  mean  pro¬ 
file  between  34°  N,  140°  W  and  32°  N,  150°  W  (refer  to  the  solid  line  in  Figure  5.8(b)  for  the 
shape  of  AC  (z) ).  This  modification  was  justified  by  the  fact  that  the  sound  speed  profiles  collect¬ 
ed  in  this  track  were  off  the  signal  propagation  path  (refer  to  Figure  2.6).  The  simulated  “acoustic 
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data”  were  generated  with  an  adiabatic  normal  mode  model  using  the  sound  speed  profiles  reflect¬ 
ing  (5.17),  while  the  replica  vectors  were  generated  using  the  original  measured  sound  speed  pro¬ 
files.  As  in  section  4.5,  we  ignored  the  bottom  effect  due  to  long-range  propagation;  only  the  first 
16  non-bottom  interacting  modes  were  used  in  the  calculation  of  the  replica  vector  and  the  simu¬ 
lated  data.  Throughout  this  chapter,  the  MUSIC  matched-field  processor  was  used  to  compute  the 
ambiguity  surfaces  due  to  its  high-resolution  capability.  For  comparison,  both  ambiguity  surfaces 
for  the  ideal  case  (no  mismatch)  and  the  mismatched  case  during  the  three  time  intervals  are  plot¬ 
ted  in  Figure  5.2  and  Figure  5.3.  respectively.  The  A’s  are  the  “true”  source  locations  according  to 
the  source  log,  and  the  *’s  are  the  MFP  estimated  source  locations  i.e.,  the  highest  peak  in  the  am¬ 
biguity  surfaces.  As  expected,  the  estimated  source  locations  using  the  mismatched  environment 
were  off  in  range  by  roughly  100  km  to  150  km  (2  to  3  CZs)  for  all  three  time  intervals. 

5.4.1. 1  Sound  Speed  Adaptation 

We  now  proceed  to  the  environment  adaptation  technique  and  assume  that  the  source-array 
geometry  is  known  exactly  during  the  first  time  interval  (01:47).  The  first  simulation  case  is  to  in¬ 
vert  for  the  range-independent  sound  speed  profile  representing  the  environment  along  the  signal 
propagation  path  between  140°  W  and  150°  W  (refer  to  Figure  2.6),  in  a  matched-field  sense.  This 
range-independent  approximation  is  justified  by  the  fact  that  this  part  of  the  ocean  is  environmen¬ 
tally  benign,  and  the  range  dependency  is  relatively  weak  as  seen  in  Figure  4.5. 

Figure  5.4  (a)  shows  the  excess  (demeaned)  sound  speed  profiles  derived  from  30  AXBT  mea¬ 
surements  made  between  34°  N,  140°  W  and  32°  N,  150°  in  the  Northeast  Pacific  in  July,  1989. 
Table  5.1  lists  the  five  largest  eigenvalues  obtained  through  eigen-decomposition  of  the  sound 

Table  5,1  The  five  largest  eigenvalues  obtained  from  eigen-decomposition  of  the  sound  speed 

covariance  matrix. 


Number 

Eigenvalue 

1 

54.06 

2 

35.96 

3 

11.93 

4 

2.81 

5 

1.81 

speed  covariance  matrix.  As  can  be  seen,  only  the  first  few  eigenvalues  are  significant.  We  thus 
use  the  eigenvectors  or  EOFs  corresponding  to  the  two  latest  eigenvalues  in  spanning  the  search 
spaces.  Normalized  versions  of  the  EOFs  corresponding  to  the  two  largest  eigenvalues  are  shown 
as  solid  and  dotted  curves,  respectively,  in  Figure  5.4(b).  The  eigenvalues  or  variances  corre- 


Figure  5.4  (a)  Excess  (demeaned)  sound  speed  profiles  computed  from  30 
AXBT  measurements  made  between  140°  W  and  150°  W,  July  1989,  (b) 
Normalized  versions  of  the  first  and  second  EOFs  derived  from  (a). 
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sponding  to  these  EOFs  were  X  j  =  o2}  =  54.06  and  X2  =  =  35.96,  respectively.  To  ensure 

all  possible  parameter  values  are  reachable  in  the  search  spaces,  we  allow  the  values  of  the  EOF 
coefficients  to  vary  within  ±60^.  Using  the  mean  sound  speed  profile,  the  first  two  EOFs,  and  the 
initial  solution  of  the  EOF  coefficients  of  (0,0),  the  optimization  is  performed  with  the  fast  simu¬ 
lated  annealing  procedure.  Figure  5.5  shows  the  convergence  properties  of  a  typical  annealing 
run.  In  this  example,  the  temperature  was  initialized  at  200  and  was  reduced  to  0.2  over  1000  iter¬ 
ations.  The  first  and  second  panels  are  the  trajectories  of  the  EOF  coefficients  while  the  third  pan¬ 
el  is  the  cost  function  learning  curve  as  the  annealing  proceeds.  As  expected,  the  cost  function 
learning  curve  declines  as  the  iteration  goes  on,  but  occasionally  the  curve  increases  to  a  higher 
energy  state  thus  indicating  the  escaping  out  of  the  local  minima.  After  350  iterations,  the  optimi¬ 
zation  process  converges  to  the  minimum  energy  state,  and  the  solution  of  the  EOF  coefficients 
found  at  the  1000th  iteration  is  (9.00, 21.91).  To  validate  the  optimal  solution,  the  energy  (cost 
function)  surface  as  a  function  of  the  EOF  coefficient  values  is  calculated  exhaustively  on  a  regu¬ 
lar  grid  in  c^and  a2(with  a  grid  size  of  1)  bounded  by  [-50,50],  which  corresponds  to  10201  eval¬ 
uations  of  the  cost  function  as  displayed  in  Figure  5.6.  The  surface  shows  a  global  minimum  at 
ctj  =  9  ,  a2  =  22  with  numerous  local  minima  scattering  around.  The  joint  trajectories  of  the 
EOF  coefficients  as  the  annealing  proceeds  are  overlaid  on  the  energy  surface  plot.  A  good  agree¬ 
ment  between  the  solutions  obtained  by  grid  search  and  optimization  is  observed.  It  is  of  interest 
to  know  whether  the  optimal  solution  of  the  EOF  coefficients  varies  from  run  to  run.  The  environ¬ 
ment  adaptation  procedure  is  then  repeated  for  nine  times,  each  time  with  1000  iterations.  The 
joint  trajectories  of  the  EOF  coefficients  for  the  nine  runs  are  plotted  in  Figure  5.7,  and  the  corre¬ 
sponding  solutions  are  given  in  Table  5.2.  The  convergency  properties  and  the  consistent  agree- 


Table  5.2  Optimization  results  from  nine  runs  using  sound  speed  EOF  coefficients  as  search 

parameters. 


Run 

EOF  Coef.  #1 

EOF  Coef.  #2 

1 

9.01 

21.88 

2 

8.96 

21.85 

3 

8.72 

22.00 

4 

8.92 

21.89 

5 

9.09 

21.93 

6 

9.05 

21.88 

7 

9.00 

21.97 

8 

9.05 

21.74 

9 

8.82 

22.02 

ment  among  all  runs  are  again  observed,  which  confirms  the  fitness  of  the  global  optimization 
algorithm. 


Normalized  Energy  .  EOPCoofW  IHFCoef#! 
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Figure  5.5 


Trajectories  of  the  sound  speed  EOF  coefficients  and  cost  function  learning 
curve  for  a  typical  annealing  run. 
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Figure  5.7  Joint  trajectories  of  sound  speed  EOF  coefficients  for  different  annealing  runs. 
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The  adapted  sound  speed  profile  derived  from  the  mean  profile  and  the  two  EOF  coefficients 
found  at  the  1000th  iteration  is  plotteu  by  dotted  line  in  Figure  5.8.  For  comparison,  the  modified 
sound  speed  profiles  (i.e.,  the  ‘’true”  profiles)  used  to  simulate  the  “acoustic  data”  are  averaged 
and  plotted  by  solid  line,  and  the  measured  profiles  that  were  used  to  model  the  replica  vectors  are 
averaged  and  plotted  by  dashed  line.  As  can  be  seen,  the  adapted  profile  tracks  the  “true”  profile 
closely.  We  now  proceed  to  the  localization  phase  of  the  technique.  The  measured  sound  speed 
profiles  between  140°  W  and  150°  W  are  replaced  by  the  single  adapted  sound  speed  profile,  and 
normal  matched-field  processing  resumes.  Figure  5.9  shows  the  optimized  range-azimuth  ambi¬ 
guity  surfaced  for  all  three  time  intervals.  The  optimized  (or  self-cohered)  source  locations  match 
those  of  an  ideal  simulation  (no  mismatch). 

5.4.1. 2  Wave  Number  Adaptation 

The  second  implementation  of  the  environmental  adaptation  technique  is  to  invert  for  the 
modal  wave  numbers  at  the  source  frequency  (14  Hz)  for  the  acoustic  waveguide  of  interest.  The 
assumptions  and  conditions  are  identical  to  the  sound  speed  adaptation  case  with  one  exception, 
that  is,  the  modal  wave  number  EOF  coefficients  are  now  the  search  parameters.  To  compare 
these  two  implementations,  we  followed  exactly  the  same  path  of  investigation.  Figure  5.10(a) 
shows  the  excess  (demeaned)  wave  number  estimates  of  the  modal  wave  number  at  the  source 
frequency  (14  Hz)  obtained  from  30  AXBT  measurements  made  in  July,  1989.  Table  5.3  lists  the 
five  largest  eigenvalues  obtained  through  eigen-decomposition  of  the  wave  number  covariance 
matrix.  Normalized  versions  of  the  EOFs  corresponding  to  the  two  largest  eigenvalues  are  shown 


Table  5.3  The  five  largest  eigenvalues  obtained  from  eigen-decomposition  of  the  wave  number 

covariance  matrix. 


Number 

Eigenvalue 

1 

1.43  X  10"V 

2 

1.06  X  10-10 

3 

4.06  X  10“" 

4 

4.55  X  10-12 

5 

2.00  X  10-’2 

as  solid  and  dashed  curves  in  Figure  5.10(b).  Note  that  the  lag  in  the  wave  number  covariance  ma¬ 
trix  estimate  is  now  the  mode  instead  of  the  depth.  Optimization  is  performed  using  the  fast  simu¬ 
lated  annealing  algorithm.  Figure  5.11  shows  the  trajectories  of  the  wave  number  EOF 
coefficients  and  the  cost  function  learning  curve  for  a  typical  run.  Figure  5.12  shows  the  energy 
surface  computed  by  exhaustive  search  and  the  joint  trajectories  of  the  wave  number  EOF  coeffi- 
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Figure  5.8  (a)  Sound  speed  profile  derived  from  the  environment  adaptation  technique 
(dotted  line),  “true”  profile  (solid  line),  and  measured  sound  speed  profile  (dashed  line); 
(b)  the  difference  version  of  (a),  i.e.,  using  the  measured  profile  as  reference  (dashed  line), 
the  “true”  minus  the  measured  and  the  adapted  minus  the  measured  are  in  solid  and  dotted 

lines,  respectively  (simulation). 
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Figure  5.10  (a)  Excess  (demeaned)  horizontal  wave  numbers  at  1 4  Hz 
derived  form  30  AXBT  measurements  made  tx  ‘ween  140°  W  and  150°  W, 
July  1989,  (b)  Normalized  versions  of  the  first  and  second  EOFs  derived 

from  (a). 
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Figure  5.12  Energy  surface  as  a  function  of  wave  number  EOF  coeffici 

by  exhaustive  search  (simulation). 
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cients  as  the  annealing  proceeds.  Figure  5.13  presents  the  trajectories  for  the  nine  runs  while  Table 
5.4  lists  the  optimal  wave  number  EOF  coefficients.  Again,  good  convergence  properties  and 


Table  5.4  Optimization  results  from  nine  runs  using  wave  number  EOF  coefficients  as  search 

parameters. 


Run 

EOF  Coef.  #1 

EOF  Coef.  #2 

1 

20.05 

-16.66 

2 

19.99 

-16.85 

3 

19.98 

-16.85 

4 

19.93 

-16.83 

5 

20.07 

-16.60 

6 

20.01 

-16.72 

7 

19.95 

-16.85 

8 

20.00 

-16.83 

9 

20.02 

-16.77 

agreement  among  the  wave  number  EOF  coefficient  estimates  are  observed.  Using  as  a  reference 
the  averaged  model  wave  numbers  derived  from  the  measured  profiles,  the  “true”  minus  the  mea¬ 
sured  and  the  adapted  minus  the  measured  are  plotted  by  solid  and  dotted  lines,  respectively,  in 
Figure  5.14.  The  adapted  wave  numbers  track  the  “true”  wave  numbers  nicely  except  for  the  first 
five  modes.  A  possible  explanation  of  the  deviation  from  the  “true”  wave  numbers  is  that  these 
modes  are  weakly  excited,  thus  less  weight  is  given  during  the  adaptation  process.  Figure  5.15 
shows  the  range-azimuth  ambiguity  surface  for  all  three  time  intervals  using  the  adapted  wave 
numbers.  The  optimized  (or  self-cohered  source)  locations  match  those  of  ideal  simulation.  It  is  of 
interest  to  know  whether  the  true  source  location  is  unique  in  the  neighborhood  of  the  source  loca¬ 
tion.  The  adaptation  procedure  is  then  repeated  for  each  assumed  source  location  within  a  spatial 
window  extended  16  km  in  range  and  1.5  degrees  in  azimuth  (approximately  60  km  in  cross 
range)  enclosing  the  true  source  location.  Figure  5.16  confirms  that  the  range-azimuth  maximum 
energy  surface  indeed  has  a  maximum  that  corresponds  to  the  true  source  location.  This  suggests 
that  a  weakly  range-dependent  environment  can  be  approximated  by  a  range-independent  envi¬ 
ronment  in  a  matched-field  sense,  given  that  the  sensor-array  geometry  is  known  exactly. 


5.4. 1. 3  Remarks 

Although  both  implementations  of  the  environment  adaptation  technique  produce  similar  re¬ 
sults  under  the  example  given  above,  the  computation  burden  for  the  two  approaches  differs  sig¬ 
nificantly.  For  the  sound  speed  case,  evaluation  of  the  cost  function  requires  the  invocation  of  the 
acoustic  propagation  model  to  output  the  modal  eigen-functions  and  wave  numbers  (this  is  com- 
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Figure  5.14  Deviations  in  adapted  wave  numbers  and  “true”  wave  numbers  from  the 
measured  model  wave  numbers  (simulation). 
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Figure  5.16  Range-azimuth  energy  surface  derived  from  repeating  the 
optimization  procedure  for  each  range-azimuth  cell  (simulation). 
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putation  intensive).  In  the  wave  number  adaptation  case,  the  wave  numbers  are  readily  available, 
and  updating  the  wave  numbers  according  to  the  perturbations  involves  only  a  few  trivial  algebra¬ 
ic  steps.  Thus,  the  difference  in  computation  performance  for  the  two  implementations  is  at  least 
on  the  order  of  a  few  magnitudes. 

5.4.2  Experimental  Data 

Data  collected  by  the  Swallow  float  array  during  the  July  1989  experiment  in  the  Northeast 
Pacific  were  then  processed  in  the  same  fashion  as  the  simulated  data.  The  acoustic  modeling  was 
the  same  as  in  the  simulation.  The  data  considered  here  were  recorded  at  01:47, 04:38,  and  07:32 
PST  July  9, 1989.  A  2048-point  data  block  for  each  time  interval  was  extracted  to  estimate  the  ar¬ 
ray  cross-spectral  density  matrix.  Each  time  interval  was  41  seconds  long.  Data  were  processed 
with  128-point  fast  Fourier  transforms  with  50  percent  overlap.  Thirty-one  snapshots  were  ob¬ 
tained  to  estimate  the  array  cross-spectral  density  matrix.  Figure  5.17  plots  the  matched-field 
range-azimuth  ambiguity  surfaces  for  all  three  intervals  with  the  highest  peak  marked  with  *’s 
and  the  true  source  locations  marked  with  A’s.  The  shift  in  the  source  locations  was  due  to  the  en¬ 
vironmental  mismatch  as  diagnosed  through  simulation  in  the  last  section.  We  then  entered  into 
the  adaptation  phase  of  the  technique.  We  cautiously  selected  the  reference  source  location  that 
corresponded  to  01:47  PST  July  9,  1989,  by  repeating  the  adaptation  procedure  for  all  assumed 
range-azimuth  cells  in  the  neighborhood  of  the  source  location  gotten  from  the  source  log.  The 
wave  number  EOF  coefficients  were  used  as  the  search  parameters  for  computational  efficiency. 
Figure  5.18  shows  the  range-azimuth  maximum  energy  surface.  The  highest  peak  in  the  surface  (a 
range  of  2489  km  and  an  azimuth  172.1°)  was  used  as  the  reference  location. 

5.4.2. 1  Sound  Speed  Adaptation 

First,  the  sound  speed  EOF  coefficients  were  used  as  the  search  parameters.  Figure  5. 19  shows 
the  trajectories  of  the  sound  speed  EOF  coefficients  during  an  annealing  run  with  the  energy  sur¬ 
face  by  grid  search  as  the  background.  Figure  5.20  plots  the  adapted  sound  speed  profile  derived 
from  the  optimal  EOF  coefficients.  The  extensive  excursion  of  the  adapted  profile  between  100  to 
200  meters  suggests  that  the  second  EOF  coefficient  plays  a  major  role  in  forming  the  underlying 
environment.  Examining  the  energy  surface  confirms  such  conjecture,  i.e.,  the  surface  peaks  at 
around  (-6, 29).  The  adapted  sound  speed  profile  was  then  used  to  replace  the  measured  sound 
speed  profiles  between  140°  W  and  150°  W.  We  then  proceeded  to  the  localization  phase  of  the 
adaptation  technique.  Figure  5.21  shows  the  optimized  range-azimuth  ambiguity  surface  for  all 
three  time  intervals.  As  can  be  seen,  the  optimized  or  self-cohered  source  track  mimics  the  track 
derived  from  the  source  log.  The  minor  discrepancy,  a  0.6°  shift  between  the  true  and  estimated 
source  locations,  is  thought  to  be  due  to  the  slight  uncertainty  in  selecting  the  coordinate  system 
with  respect  to  true  north  for  the  float  localization  (refer  to  Figure  3.14),  since  the  orientation  of 
the  X  axis  of  the  coordinate  system  was  taken  to  be  the  ship’s  position  when  bottomed  floats  9  and 
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Figure  5.17  Matched-field  processing  of  experimental  data  during  the 

three  time  intervals. 
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Figure  5.18  Range-azimuth  energy  surface  derived  from  repeating  the 
optimization  procedure  for  each  range-azimuth  cell  using  data  collected  at 

01:47  PST,  July  1989. 
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Figure  5.20  (a)  Sound  speed  profile  (dotted  line)  computed  from  the  optimal  EOF 
coefficients  obtained  from  Figure  5.19  and  the  averaged  model  sound  speed  profile 
(dashed  line),  (b)  The  difference  version  of  (a)  i.e.,  using  the  measured  profile  as 
reference  (dashed  line),  the  adapted  minus  the  measured  is  plotted  in  dotted  line. 
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10  were  put  into  the  water.  A  relative  motion  of  approximating  60  meters  in  the  Y  direction  be¬ 
tween  floats  9  and  10  while  the  floats  descended  to  the  bottom  would  cause  the  0.6°  rotation  in  the 
MFP  results.  This  order  of  error  between  ship  position  and  true  float  position  on  the  bottom  has 
been  noted  in  other  Swallow  float  experiments. 

S.4.2.2  Wave  Number  Adaptation 

We  next  used  the  wave  number  EOF  coefficients  as  the  search  parameters.  Figure  5.22  shows 
the  trajectories  of  the  wave  number  EOF  coefficients,  corresponding  to  the  annealing  run  that  pro¬ 
duced  the  highest  peak  in  Figure  5.18.  The  energy  surface  is  plotted  as  the  background.  Figure 
5.23  plots  the  first  16  modal  -iave  numbers  (derived  from  the  wave  number  EOF  coefficients) 
with  respect  to  the  measured  wave  numbers.  The  adapted  modal  wave  numbers  were  used  to  re¬ 
place  the  measured  modal  wave  numbers  between  140°  W  and  150°  W.  We  then  proceeded  to 
the  localization  phase  of  the  adaptation  technique.  Figure  5.24  shows  the  optimized  range-azi¬ 
muth  ambiguity  surface  for  all  three  time  intervals  for  the  wave  number  adaptation  approach.  The 
optimized  (self-cohered)  source  locations  match  those  of  the  sound  speed  approach  except  that 
the  source  location  during  the  third  interval,  i.e.,  07:32,  is  missed  by  one  CZ.  The  range  estimate 
error  found  during  the  third  time  interval  is  thought  to  be  due  to  the  fact  that  the  wave  number  ad¬ 
aptation  approach  adapts  the  wave  numbers  only  and  is  not  robust  enough  to  accommodate  large 
variations  in  the  modal  depth  eigenfunctions  at  the  source,  i.e.,  the  assumed  modal  eigenfunctions 
at  the  source  differ  significantly  from  the  true  underlying  modal  eigenfunctions  at  the  source. 

5.5  Summary 

The  major  difficulty  with  MFP  is  the  environmental  mismatch  problem.  In  this  chapter,  two 
techniques  to  deal  with  this  problem  were  reviewed.  An  approach  to  the  environmental  mismatch 
problem  was  presented.  The  key  idea  of  this  approach  was  the  adaptation  of  the  replica  pressure 
fields  to  the  acoustic  data  received  by  the  sensor  array  so  that  the  environmental  parameters  along 
the  path  between  a  reference  source  and  the  sensor  array  can  be  inverted  assuming  the  source-ar¬ 
ray  geometry  is  known  exactly.  Using  the  inverted  (or  adapted)  environmental  parameters,  normal 
MFP  can  be  resumed  to  search  for  an  unknown  source  of  interest.  An  efficient  optimization  algo¬ 
rithm  was  implemented  using  the  EOF  approach  to  compress  the  search  spaces  and  a  fast  anneal¬ 
ing  algorithm  to  search  the  optimal  environmental  parameter  values.  Two  implementations  of  the 
environmental  adaptation  technique  were  first  studied  through  simulation  and  then  applied  to  ex¬ 
perimental  data.  Both  simulation  and  real  data  processing  results  demonstrated  the  potential  of 
this  technique.  The  self-cohered  14-Hz  source  locations  for  a  period  of  6  hours  during  the  July 
1989  experiment  were  found  to  be  consistent  with  the  source  locations  obtained  from  the  source 
log.  The  sound  speed  structure  adaptation  implementation  of  the  technique  seemed  to  be  robust 
but  computation  cumbersome.  However,  the  modal  wave  number  adaptation  implementation  was 
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168 


0.0 


26 


computationally  efficient  but  sensitive  to  separation  between  the  reference  and  the  unknown 
source  locations  in  time  and  space. 
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CONCLUSIONS 


The  central  theme  of  this  dissertation  has  been  the  underwater  source  localization  prob¬ 
lem.  Specifically,  we  have  demonstrated  the  source  localization  and  tracking  capability  of  a  freely 
drifting  volumetric  array  with  MFP  using  experimental  data.  We  have  also  proposed  and  demon¬ 
strated  an  environment  adaptation  technique  to  deal  with  the  environmental  mismatch  problem. 

Data  collected  during  the  19S9  Swallow  float  experiment  were  used  to  perform  the  study. 
The  geometries  of  the  Swallow  float  array  as  a  function  of  time  during  the  1989  experiment  have 
been  estimated  using  the  8-kHz  range  measurement  with  both  least-squares  and  Kalman  filter 
methods.  The  rms  position  errors  estimated  by  both  methods  are  within  5  meters,  which  is  within 
the  desired  accuracy  of  one-tenth  of  a  wavelength  at  the  highest  frequency  of  interest  25  Hz  (6  m) 
in  order  to  effectively  process  the  VLF  acoustic  data  coherently.  Furthermore,  analysis  of  the  ex¬ 
periment  acoustic  data  showed  high  signal-to-noise  ratio  and  high  coherence  at  14  Hz  among  all 
nine  freely  drifting  floats  during  some  time  intervals.  The  14-Hz  was  a  continuous  wave  tonal  pro¬ 
jected  by  the  Mark  VI  from  the  R/V  Aloha  involving  in  a  companion  experiment.  The  high  coher¬ 
ence  among  the  floats  provided  an  opportunity  for  matched-field  beamforming  of  the  VLF 
acoustic  data. 

The  replica  vectors  were  modeled  with  an  adiabatic  normal  mode  numerical  technique  us¬ 
ing  the  environmental  data  collected  during  the  experiment.  Initial  matched-field  processing  of 
the  experimental  data  experienced  difficulties  in  estimating  the  source  depth  and  range  while  the 
source  azimuth  estimate  was  somewhat  successful.  Controlled  simulation  using  the  same  condi¬ 
tions  as  in  the  real  data  has  revealed  that  (1)  depth  resolution  indeed  is  a  difficult  problem  for  a 
shallow  VLF  source  in  a  long-range  environment  due  to  fewer  modes  being  available  and  due  to 
low  source  depth-to-wavelength  ratio,  (2)  the  range  estimate  is  sensitive  to  environmental  mis¬ 
match,  and  (3)  the  azimuth  estimate  is  robust.  The  main  cause  of  the  performance  degradation  has 
thus  been  identified  to  be  uncertainty  in  the  environment  (i.e.,  sound  speed  mismatch). 

An  environment  adaptation  technique  using  a  global  optimization  algorithm  was  devel¬ 
oped  to  counteract  the  MFP  performance  degradation  due  to  uncertainty  in  the  ocean  acoustic  en¬ 
vironment.  We  have  demonstrated  through  simulation  that  with  limited  a  priori  knowledge  of  the 
environment  and  with  a  reference  source  at  a  known  location,  the  environment  can  be  adapted  in  a 
matched-field  sense.  Using  the  adapted  environment,  other  unknown  source(s)  of  interest  in  the 
vicinity  of  the  reference  source  can  be  correctly  localized.  Applying  the  environment  adaptation 
technique  to  experimental  data  has  shown  that  the  14-Hz  source  was  successfully  localized  and 
tracked  in  azimuth  and  range  within  a  region  of  interest  using  the  MFP  technique  for  a  period  of  6 
hours. 
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There  are  several  directions  in  which  the  results  developed  here  can  be  extended.  In  this 
dissertation,  due  to  the  nature  of  the  data  set,  emphasis  was  placed  on  localizing  underwater  sound 
source  at  long  range.  It  would  be  of  interest  to  further  test  the  localization  performance  of  the 
Swallow  float  array  and  the  environment  adaptation  concept  in  intermediate  and  close -range  situ¬ 
ations  using  experimental  data.  Also,  we  need  to  quantify  the  performance  limitations  of  the  envi¬ 
ronment  adaptation  technique  through  extensive  simulation. 

This  thesis  has  been  restricted  to  using  the  acoustic  pressure  data  from  the  hydrophones. 
Since  acoustic  particle  velocity  data  from  three  channels  of  the  geophones  are  also  available,  it 
would  be  of  interest  to  process  the  particle-velocity  data  in  a  matched-field  sense  and  compare  the 
results  to  those  of  acoustic-pressure  data. 

Finally,  although  the  environment  adaptation  technique  did  perform  well,  efforts  should 
also  be  directed  toward  investigating  a  more  robust  array  processing  technique,  i.e.,  developing  an 
MFP  processor  that  is  tolerant  to  some  degrees  of  mismatch  between  replica  and  data  and  that  still 
maintains  a  satisfactory  localization  performance. 
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APPENDICES 


Harmonic  Mean  Sound  Speed  Estimates 


Floats 

Sound  Speed 
(m/msec) 

Error  Variance 

(  (m/msec  ) 

^  9  &  10 

1.5237 

0.1  X  lO"-6 

9  &  11 

1.5237 

0.1  X  10"6 

9  &  0 

1.4975 

0.1  X  lO-* 

9  &  1 

1.4999 

0.1  X  10-6 

9  &  2 

1.5024 

0.1  X  10-6 

9  &  3 

1.5043 

0.1  X  10-6 

9  &  4 

1.5075 

0.1  X  10~* 

9  &  5 

1.5108 

0.1  X  10-6 

9  &  6 

1.5143 

0.1  X  10-6 

9  &  7 

1.5178 

0.1  X  10~* 

9  &  8 

1.5213 

o.i  x  io-6 

10&  11 

1.5237 

0.1  X  10-6 

10&0 

1.4975 

0.1  X  10-6 

10  &  1 

1.4999 

0.1  X  lO"6 

10  &  2 

1.5024 

0.1  X  10-* 

10  &  3 

1.5043 

0.1  X  10~* 

10  &  4 

1.5075 

0.1  X  10-6 

10  &  5 

1.5108 

0.1  X  10-6 

10  &  6 

1.5143 

O.I  X  10~* 

10&  7 

1.5178 

0.1  X  lO-6 

10  &  8 

1.5213 

0.1  X  10-6 

11  &0 

1.4975 

0.1  X  lO--6 

11  &  1 

1.4999 

0.1  X  10-6 

11  &  2 

1.5024 

0.1  X  10"6 

11  &  3 

1.5043 

0.1  X  lO"6 

11  &  4 

1.5075 

0.1  X  10""6 

11  &  5 

1.5108 

0.1  X  10"6 

11  &  6 

1.5143 

0.1  X  10-6 

11  &  7 

1.5178 

0.1  X  10-6 

11  &  8 

1.5213 

0.1  X  10-6 

0&  1 

1.4809 

0.1  X  10-6 

0  &  2 

1.4820 

0.1  X  10-6 

0  &  3 

1.4830 

o.i  x  io-6 

0  &  4 

1.4849 

0.1  X  10~6 

0  &  5 

1.4872 

0.1  X  10"-6 

0  &  6 

1.4897 

O.I  X  10-6 

0  &  7 

1.4926 

0.1  X  10~6 

Floats 

Sound  Speed 
(m/ msec) 

Error  Variance 

(  (m/msec  )  “  ) 

0  &  8 

1.4955 

0.1  X  10"* 

1  &2 

1.4832 

0.1  X  10~* 

1  &3 

1.4844 

0.1  X  10~* 

1  &  4 

1.4865 

0.1  X  10"* 

1  &  5 

1.4890 

0.1  X  10~* 

1  &6 

1.4917 

0.1  X  10"* 

1  &  7 

1.4948 

0.1  X  10~* 

1  &  8 

1.4978 

0.1  X  10~* 

2  &  3 

1.4860 

0.1  X  10"* 

2  &  4 

1.4883 

0.1  X  10"* 

2  &  5 

1.4909 

0.1  X  10"* 

2  &  6 

1.4939 

0.1  X  10“* 

2  &  7 

1.4970 

0.1  X  !0"* 

2  &  8 

1.5002 

0.1  X  I0~* 

3  &  4 

1.4898 

0.1  X  10"* 

3  &  5 

1.4926 

0.1  X  10"* 

3  &  6 

1.4957 

0.1  X  10"* 

3  &  7 

1.4989 

0.1  X  10"* 

3  &  8 

1.5022 

0.1  X  10"* 

4  &  5 

1.4954 

0.1  X  10"* 

4  &  6 

1.4986 

0.1  x  to"* 

4  &  7 

1.5020 

0.1  X  10-6 

4  &  8 

1.5053 

0.1  X  10"* 

5  &  6 

1.5018 

0.1  X  10"* 

5  &  7 

1.5052 

0.1  X  10"* 

5  &  8 

1.5086 

0.1  X  10"* 

6  &  7 

1.5086 

0.1  X  10"* 

6  &  8 

1.5120 

0.1  X  10"* 

7  &  8 

1.5155 

0.1  X  10"* 

9  &  9 

0.7478 

0.1  X  10"* 

10&  10 

0.7478 

0.1  X  10"* 

11  &  11 

0.7478 

0.1  X  10"* 

O&O 

0.7423 

0.1  X  10"* 

1  &  1 

0.7415 

0.)  X  10"* 

2  &  2 

0.7415 

0.1  X  10"* 

3  &  3 

0.7418 

0.1  X  10"* 

4  &  4 

0.7424 

0.1  X  10"* 

5  &  5 

0.7433 

0.1  X  10"* 

6  &  6 

0.7443 

0.1  X  10"* 

7  &  7 

0.7456 

0.1  X  10"* 

_ 

0.7469 

0.1  X  10"* 
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B:  Reciprocal  Path  Travel  Time  Variances 


Floats 

Variance 

(msec2) 

1  &  6 

2.4348 

1  &7 

1.5548 

1  &  8 

1.0353 

2  &  3 

0.9688 

2  &  4 

0.8213 

2  &  5 

1.5815 

2  &  6 

1.4629 

2  &  7 

1.3690 

2  &  8 

1.5648 

3  &  4 

0.3759 

3  &  5 

0.7077 

3  &  6 

2.2006 

3  &  7 

0.8939 

3  &  8 

2.3961 

4  &  5 

0.4745 

4  &  6 

0.3362 

4  &  7 

1.1722 

4  &  8 

2.3950 

5  &  6 

0.4461 

5  &  7 

0.5171 

5  &  8 

0.2398 

6  &  7 

0.9575 

6  &  8 

0.7759 

7  &  8 

0.8713 
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C:  Coefficients  of  Second-Order  Fits  to  the  Reciprocal  Path  Measurement 
Differences 


Float 

Float 

Curvature 

(msec/(100 

records)2) 

Slope  (msec/ 
100  records) 

Intercept 

(msec) 

9 

10 

0.01518 

-2.50139 

-44.85477 

9 

11 

0 

54.54246 

-70.34302 

9 

0 

-0.00399 

25.55177 

-143.75967 

9 

1 

-0.02566 

36.63271 

-180.84509 

9 

2 

0.10133 

-0.61512 

-120.54930 

9 

3 

0.07983 

42.33369 

-39.60468 

9 

4 

0.07021 

34.60133 

4.20795 

9 

5 

0.10593 

24.74624 

-73.47556 

9 

6 

0.00936 

10.18588 

-74.90748 

9 

7 

-0.00872 

-9.31168 

-56.46815 

9 

8 

-0.05318 

21.87121 

-23.24966 

10 

11 

-0.04577 

58.26873 

-36.45166 

10 

0 

0 

27.24289 

-89.77655 

10 

1 

0.00429 

38.10869 

-131.25516 

10 

2 

0.01703 

3.83972 

-86.58073 

10 

3 

0.04859 

45.28195 

3.27728 

10 

4 

0 

38.47172 

40.12939 

10 

5 

0 

29.46763 

-41.03638 

10 

6 

-0.04374 

13.62054 

-33.69858 

10 

7 

0 

-7.32420 

-10.27667 

10 

8 

0.00109 

22.64396 

30.63046 

11 

0 

0 

-29.19710 

-73.01617 

11 

1 

0.05370 

-19.95675 

-101.24307 

11 

2 

0.05198 

-53.80684 

-58.25610 

11 

3 

0.09517 

-12.82133 

35.28806 

11 

4 

0.04854 

-19.51229 

69.53906 

11 

5 

0.18019 

-31.36666 

-0.52029 

11 

6 

0.01447 

-44.42154 

-6.89490 

11 

7 

0.02574 

-64.76017 

16.32593 

11 

8 

0.01572 

-35.06775 

64.42834 

0 

1 

0.02894 

9.89287 

-32.20901 

0 

2 

0.02596 

-23.89109 

8.72311 

0 

3 

0.01238 

18.85915 

91.22421 

0 

4 

0.04712 

10.02173 

139.38211 

0 

5 

0.08044 

0.14604 

60.69196 

0 

6 

-0.04036 

-13,76737 

57.26703 

0 

7 

0.01477 

-35.07599 

84.37808 
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Curvature 

(msec/(100 

records)2) 

Slope  (msec/ 
100  records) 

Intercept 

(msec) 

0.00207 

-4.94819 

126.70593 

0.02212 

-34.60089 

47.02808 

-0.04047 

9.60653 

118.90907 

-0.02242 

1.17007 

165.56229 

0.05049 

-9.63175 

92.00610 

-0.08490 

-23.35032 

88.31992 

-0.03788 

-44.52425 

114.87286 

-0.03449 

-14.65317 

158.19537 

0.01816 

41.85386 

86.79614 

-0.12803 

37.56519 

109.27396 

0.02849 

24.51210 

50.82257 

-0.01694 

8.64918 

58.06210 

0.00913 

-11.82151 

79.21877 

-0.09417 

20.80200 

107.27597 

-0.04826 

-6.73761 

36.29120 

0.08602 

-19.59419 

-22.  :C524 

-0.04837 

-32.81397 

-32.02466 

-0.06930 

-52.08339 

-16.81482 

-0.02307 

-23.57817 

35.45761 

0.07278 

-10.60359 

-74.41269 

0 

-26.25876 

-64.9363’ 

0 

-45.90385 

-51.37494 

-0.06807 

-14.64013 

-13.09142 

-0.12665 

-13.52716 

-5.83221 

-0.05768 

-35.12783 

20.22412 

-0.05973 

-5.68480 

69.97832 

0.05221 

-21.31087 

26.96481 

-0.03280 

10.73389 

59.33090 

0.07448 

27.69873 

58.11115 
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D:  Magnitude-Squared  Coherence  (Record  1145) 


Float 

MSC 

1 

0.750720 

2 

0.885170 

3 

0.700810 

4 

0.712970 

5 

0.804030 

6 

0.679740 

7 

0.820220 

8 

0.816530 

2 

0.789920 

3 

0.616370 

4 

0.580990 

5 

0.634840 

6 

0.627240 

7 

0.666920 

8 

0.670200 

3 

0.712980 

4 

0.755040 

5 

0.774810 

6 

0.802950 

7 

0.814410 

8 

0.836930 

4 

0.609960 

5 

0.726390 

6 

0.739470 

7 

0.710770 

8 

0.742590 

5 

0.732270 

6 

0.662420 

7 

0.617800 

8 

0.683280 

6 

0.695380 

7 

0.689330 

8 

0.765470 

7 

0.638040 

8 

0.713870 

8 

0.777410 
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